REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  NO.  0704-0188 


The  public  reporting  burden  for  this  coilection  of  information  is  estimated  to  average  1  hour  per  response,  inciuding  the  time  for  reviewing  instructions, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  compieting  and  reviewing  the  coiiection  of  information.  Send  comments 
regarding  this  burden  estimate  or  any  other  aspect  of  this  coilection  of  information,  including  suggesstions  for  reducing  this  burden,  to  Washington 
Headquarters  Services,  Directorate  for  information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Ariington  VA,  22202-4302. 
Respondents  shouid  be  aware  that  notwithstanding  any  other  provision  of  iaw,  no  person  shaii  be  subject  to  any  oenaity  for  failing  to  comply  with  a  coiiection 
of  information  if  it  does  not  dispiay  a  currentiy  vaiid  OMB  controi  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


5c.  PROGRAM  ELEMENT  NUMBER 
611102 


5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 


12.  DISTRIBUTION  AVAILIBILITY  STATEMENT 
Approved  for  public  release;  distribution  is  unlimited. 

13.  SUPPLEMENTARY  NOTES 

The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and  should  not  contrued  as  an  official  Department 
of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other  documentation. 

14.  ABSTRACT 


15.  SUBJECT  TERMS 


17.  LIMITATION  OF  15.  NUMBER 
ABSTRACT  OF  PAGES 

UU 

Standard  Form  298  (Rev  8/98) 
Prescribed  by  ANSI  Std.  Z39. 18 


19a.  NAME  OF  RESPONSIBLE  PERSON 

Qingda  Yang _ 

19b.  TELEPHONE  NUMBER 
305-284-3221 


16.  SECURITY  CLASSIFICATION  OF: 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

UU 

UU 

UU 

7.  PERFORMING  ORGANIZATION  NAMES  AND  ADDRESSES 

University  of  Miami  -  Coral  Gables 
1320  South  Dixie  Highway 
Suite  650,  EC  2960 

Coral  Gables,  FL _ 33146  -2926 _ 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS 
(ES) 

U.S.  Army  Research  Office 
P.O.  Box  12211 

Research  Triangle  Park,  NC  27709-2211 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 
ARO 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

63154-EG.28 


5d.  PROJECT  NUMBER 


3.  DATES  COVERED  (From  -  To) 
15-Jul-2013  -  14-Jul-2017 

5a.  CONTRACT  NUMBER 

W911NF-13-1-0211 _ 

5b.  GRANT  NUMBER 


2.  REPORT  TYPE 
Final  Report 

4.  TITLE  AND  SUBTITLE 

Final  Report:  ARO  1.2:  Solid  Mechanics:  Augmented  Finite 
Element  Method  For  High-Fidelity  Analysis  of  Structural 
Composites 

6.  AUTHORS 


1.  REPORT  DATE  (DD-MM-YYYY) 
03-10-2017 


RPPR  Final  Report 

as  of  16-Oct-2017 


Agency  Code: 

Proposal  Number:  63154EG  Agreement  Number:  W911NF-1 3-1 -0211 

INVESTIGATOR(S): 

Name:  Ph.D  Qingda  Yang 
Email:  qdyang@miami.edu 

Phone  Number:  3052843221 
Principal:  Y 

Organization:  University  of  Miami  -  Coral  Gables 

Address:  1320  South  Dixie  Highway,  Coral  Gables,  FL  331462926 
Country:  USA 

DUNS  Number:  625174149  EIN:  590624458 

Report  Date:  14-Oct-2017  Date  Received:  03-0ct-2017 

Final  Report  for  Period  Beginning  15-Jul-2013  and  Ending  14-Jul-2017 

Title:  ARO  1.2:  Solid  Mechanics:  Augmented  Finite  Element  Method  For  High-Fidelity  Analysis  of  Structural 
Composites 

Begin  Performance  Period:  15-Jul-2013  End  Performance  Period:  14-Jul-2017 

Report  Term:  O-Other 

Submitted  By:  Ph.D  Qingda  Yang  Email:  qdyang@miami.edu 

Phone:  (305)284-3221 

Distribution  Statement:  1-Approved  for  public  release;  distribution  is  unlimited. 

STEM  Degrees:  3  STEM  Participants:  6 

Major  Goals:  The  proposed  objective  of  this  study,  according  to  the  original  propose,  is  to  develop  “an  efficient 
and  accurate  multiscale  computational  methodology  for  rapid  virtual  testing  of  composites  with  complex 
microstructures”.  The  proposed  research  goals  are: 

“1)  Development  of  new  2D  and  3D  A-FEM  elements  that  can  explicitly  account  for  statistical  material 
heterogeneity  and  can  faithfully  predict  progressive  damage  initiation,  propagation  and  multi-crack  interaction”; 

“2)  Apply  the  new  method  to  study  composites  of  interest  to  the  Army,  coupled  with  3D  micro-CT-based  precision 
experiments,  to  understand  key  physical  phenomena  of  crack  formation  and  growth,  and  to  quantify  their  direct 
effects  on  structural  integrity  under  general  loads”; 

Two  more  goals  were  added  in  the  4th  year  extension  proposal  in  view  of  the  importance  of  material  nonlinearity 
due  to  finer  scale  (micron-  or  sub-micron  scale)  material  damage, 

“3)  Complete  the  development  of  a  ‘multi-physics’  nonlinear  element  within  the  AFEM  framework  that  will  allow  us 
to  represent  a  microcrack  in  a  heterogeneous  medium”; 

“4)  Verify  the  ability  of  the  multi-physics  element  to  make  correct  predictions  of  the  conditions  for  crack  deflection, 
bifurcation,  and  branching  and  of  the  influence  on  these  phenomena  of  fine-scale  continuum  damage  in  the 
surrounding  material”. 

Accomplishments:  Major  Technical  Achievements  resulted  from  the  projects: 

•  A  new  augmentation  formulation  that  enables  standard  finite  element  method  to  deal  with  strong  and  weak 
discontinuities  caused  by  arbitrary,  coupled  multiple  damage  formation  and  propagation  in  complex  heterogeneous, 
without  the  need  for  adding  extra  DoFs  into  the  global  problem. 

•  A  novel  elemental  condensation  procedure  (consistency-check  based  algorithm)  that  demonstrated  2-3  orders  of 
magnitude  improvement  in  numerical  efficiency,  accuracy,  and  robustness. 

•  A  novel  inertia-based  stabilizing  method  that  guarantees  unconditional  convergence  and  allows  for  smooth 
transition  from  quasi-static  crack  development  to  rapid  or  even  dynamic  crack  propagation,  and  vice  versa. 

•  A  nonlinear  A-FEM  to  account  for  both  finer  scale  (i.e.  submicron  scale)  continuum  damage  or  plasticity  and 
larger  scale  nonlinearity  due  to  nonlinear  coupling  of  multiple  cracks. 

•  Extended  the  new  A-FEM  with  coupled  thermal-mechanical  analysis  capabilities  under  static,  cyclic,  and  dynamic 
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loading  conditions. 

•  Much  improved  understanding  of  the  interaction  of  multiple  damage  modes  in  composites  and  their  convoluted 
influence  on  structural  integrity  through  comparing/correlating  detailed  A-FEM  analysis  results  against  high 
resolution  experimental  measurements  /characterizations  (DIG,  micro-CT,  etc) 

Training  Opportunities:  This  project  helped  to  train  the  following  graduate  students  and  post-doctorate 
researchers: 

1.  Dr.  Mehdi  Naderi  (former  post-doc  researcher,  12/2012-4/2014),  80%  Effort.  Now  with  Technical  Data  Analysis, 
Inc,  Washington,  DC. 

2.  Dr.  Jaedal  Jung  (former  post-doc  Researcher,  7/2014-4/2017),  overall  30%  effort  since  7/2014.  Now  with 
Teledyne  Scientific  Company,  Thousand  Oaks,  CA. 

3.  Dr.  Derek  Schesser  (Ph.D  student,  graduated  in  5/2015),  50%  support  from  this  project  during  2013-2015.  Now 
a  DARPA  Science  Engineering  Technology  Advisor  (SETA). 

4.  Ms.  Brigitte  Morales  (PhD  student,  6/2014-5/2015),  50%  support.  Left  this  program  and  joined  another  group  in 
the  department  in  6/2015  after  being  awarded  with  an  NSF  scholarship. 

5.  Ms.  Yunwei  Xu  (PhD  student  in  progress),  supported  by  this  project  100%  from  8/15/2014-  8/14/2016. 

6.  Mr.  Yuhang  Yang  (PhD  student  in  progress)  supported  100%  from  8/15/2016-8/14/2017. 

7.  Mr.  Scott  Floerke  (MS  student  graduated  5/2017).  Performed  C4  blast  simulations  on  aluminum  panels  under 
joined  advising  of  the  PI  and  Dr.  Chian-fong  Yen  in  ARL. 


8.  Mr.  James  Robertson  (MS  student  graduated  5/2017).  Programmed  a  fiber/matrix  model  generator  using  python 
script  language  which  enabled  high-fidelity  RVE  analyses. 
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Results  Dissemination:  This  projected  resulted  17  published  journal  papers  and  1  book  chapter.  Another  five 
journal  publications  are  in  preparation.  Below  is  the  list: 

Published/accepted  (18): 

1.  Q.  D.  Yang  and  B.  C.  Do,  “Predicting  Damage  Evolution  in  Composites  with  Explicit  Representation  of  Discrete 
Damage  Modes”,  Handbook  of  Damage  Mechanics:  Nano  to  Macro  Scale  for  Materials  and  Structures,  editor 
George  Voyiadjis.  SpringerReference,  New  York,  2014  (book  chapter) 

2.  W.  Liu,  Q.D.  Yang*,  S.  Mohammadizadeh,  and  X.Y.  Su,  “An  Efficient  Augmented  Finite  Element  (A-FEM)  for 
Arbitrary  Cracking  and  Crack  Interaction  in  Solids”,  International  Journal  of  Numerical  Methods  in  Engineering,  99 
(438-468),  2014. 

3.  D.S.  Ling*,  L.F.  Bu,  F.B.  Tu,  Q.D.  Yang,  and  Y.M.  Chen,  “A  finite  element  method  with  mesh-separation-based 
approximation  technique  and  its  application  in  modeling  crack  propagation  with  adaptive  mesh  refinement”. 
International  Journal  of  Numerical  Methods  in  Engineering,  99  (487-521),  2014 

4.  F.B.  Tu,  D.S.  Ling*,  L.F.  Bu,  and  Q.D.  Yang,  Generalized  bridging  domain  method  for  coupling  finite  elements 
with  discrete  elements.  Computer  Methods  in  Applied  Mechanics  and  Engineering,  276  (509-533),  2014. 

5.  Q.  D.  Yang*,  D.  Schesser,  M.  Niess,  P.  Wright,  M.  Mavrogordato,  I.  Sinclair,  S.  M.  Spearing,  and  B.  N.  Cox,  “Qn 
crack  initiation  in  notched,  cross-plied  polymer  matrix  composites”.  Journal  of  the  Mechanics  and  Physics  of  Solids, 
78  (314-332).  2014. 

6.  C.  X.  Zhang,  J.  Z.  Song,  Q.  D.  Yang,  “Periodic  buckling  patterns  of  graphene/hexagonal  boron  nitride 
heterostructure”.  Nanotechnology,  25  (445401-10),  2014 

7.  Q.  Xu,  W.M.  Tao,  S.X.  Qu*,  Q.  D.  Yang,  “A  cohesive  zone  model  for  the  elevated  temperature  interfacial 
debonding  and  frictional  sliding  behavior”.  Composite  Science  and  Technology,  110  (45-52),  2015 

8.  W.  Liu,  D.  Schesser,  Q.  D.  Yang  and  D.S.  Ling,  “A  Consistency-Check  Based  Algorithm  for  Element 
Condensation  in  Augmented  Finite  Element  Methods  for  Fracture  Analysis”,  Engineering  Fracture  Mechanics,  139 
(78-97),  2015. 

9.  Y.  C.  Gu,  J.  Jung,  Q.  D.  Yang,  and  W.  Q.  Chen,  “A  New  Stabilizing  Method  for  Numerical  Analyses  with  Severe 
Local  and  Global  Instability”,  ASME  Journal  of  Applied  Mechanics,  82  (101010-1,  -12),  2015 

10.  J.  Jung,  B.  C.  Do,  and  Q.  D.  Yang,  “A-FEM  for  Arbitrary  Cracking  and  Crack  Interaction  in  Solids  under 
Thermo-Mechanical  Loadings”,  Philosophical  Transactions  A.  374:  20150282.2016 

11.  G.  Borstnar,  M.  N.  Mavrogordato,  Q.  D.  Yang,  I.  Sinclair,  and  S.  M.  Spearing,  “Crack  path  simulation  in  a 
particle-toughened  interlayer  within  a  polymer  composite  laminate”.  Composite  Science  &  Technology  133  (89-96). 
2016 

12.  A.  Vanaerschot,  B.  N.  Cox,  S.  V.  Lomov,  D.  Vandepitte,  “Multi-scale  modelling  strategy  for  textile  composites 
based  on  stochastic  reinforcement  geometry,”  Comput.  Methods  Appl.  Mech.  Engrg.  310  (906-934).  2016 

13.  M.  Naderi,  J.  Jung,  and  Q.  D.  Yang,  “A  three-dimensional  Augmented  Finite  Element  for  Modeling  arbitrary 
cracking  in  solids”.  International  Journal  of  Fracture,  192  (147-168).  2016 

14.  J.  Jung  and  Q.  D.  Yang,  “A  Two-Dimensional  Augmented  Finite  Element  for  Dynamic  Crack  Initiation  and 
Propagation”.  International  Journal  of  Fracture.  203  (41-61).  2017 

15.  G.  D.  Nian,  Y.  J.  Shan,  S.X.  Qu,  and  Q.  D.  Yang,  “Failure  analysis  of  syntactic  foams:  a  computational  model 
with  cohesive  law  and  XFEM”,  Composites  Part  B.  89  (18-26),  2016 

16.  B.  N.  Legarth  and  Q.  D.  Yang,  “Micro-mechanical  analyses  of  debonding  and  matrix  cracking  in  dual-phase 
materials”,  ASME  Journal  of  Applied  Mechanics,  83  (051006-1,  -9),  2016 

17.  J.  S.  Hu,  K.  F.  Zhang,  Q.D.  Yang,  H.  Cheng,  Y.  Yang,  An  experimental  study  on  bearing  response  of  single-lap 
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bolted  composite  interference-fit  joints.  In  press,  Composite:  Part  B.  2017 

18.  J.S.  Hu,  K.  F.  Zhang,  Q.  D.  Yang,  C.  Hui,  S.  N.  Liu,  Y.  Yang,  Fretting  behaviors  of  interface  between  CFRP  and 
coated  titanium  alloy  in  composite  interference-fit  joints  under  service  condition.  In  press.  Materials  &  Design,  2017. 

Manuscripts  in  Preparation: 

19.  Y.W.  Xu,  J.  Jung,  S.  Nojavan,  and  Q.D.  Yang,  “An  Orthotropic  A-FEM  for  Progressive  Damage  Analyses  of 
Laminar  Composites  with  Explicit  Presentation  of  All  Major  Damage  Modes,”  in  preparation  for  possible  submission 
to  Composite  A,  2017. 

20.  J.  Quek  and  B.  N.  Cox,  “Graph  theory  analysis  of  rich  fiber-scale  data  opens  a  new  route  to  high-fidelity 
simulations  of  damage  evolution  in  composites,”  in  preparation  for  submission  to  JMPS,  2017. 

21.  Y.W.  Xu,  J.  Jung,  and  Q.  D.  Yang,  “Effects  of  fine  scale  nonlinearity  on  coupled  delamination  and  ply  cracking 
in  laminated  composites,”  in  preparation  for  possible  submission  to  Composites  &  Structures,  2017. 

22.  Q.D.  Yang,  T.  A.  Parthasarathy,  O.  Sudre,  B.N.  Cox,  “Coupled  Binary  Model  (CDM)  and  Continuum  Damage 
Model  (CDM)  for  Failure  Analysis  of  Textile  CMC  Structures  Subject  to  Environmental  Degradation,  ”  in 
presparation  for  possible  submission  to  the  Journal  of  American  Ceramics  Society.  2017. 

23.  Q.D.  Yang,  Y.  H.  Yang,  Q.  Sudre,  B.N.  Cox,  “Coupled  Binary  Model  (CDM)  and  Continuum  Damage  Model 
(CDM)  for  Progressive  Failure  Analysis  of  Large  Textile  CMC  Structures,  ”  in  presparation  for  possible  submission 
to  the  Journal  of  American  Ceramics  Society.  2017. 


Honors  and  Awards:  Nothing  to  Report 
Protocol  Activity  Status: 

Technology  Transfer:  •  Key  enabling  capabilities  of  the  A-FEM  have  been  enriched  and  applied  to  several 
other  governmental  and  industrial  programs. 

•  Broad  collaborations  with  national/international  universities,  research  institutions,  and  industrial  companies. 
Details  in  Appendix  C  of  the  final  report 

PARTICIPANTS: 

Participant  Type:  PD/PI 
Participant:  Qingda  Yang 

Person  Months  Worked:  8.00  Funding  Support: 

Project  Contribution: 

International  Collaboration: 

International  Travel: 

National  Academy  Member:  N 
Qther  Collaborators: 


Participant  Type:  Consultant 
Participant:  Brian  Cox 

Person  Months  Worked:  4.00  Funding  Support: 

Project  Contnbutjon: 

Internatjonal  Collaboratjon: 

Internatjonal  Travel: 

National  Academy  Member:  N 
Qther  Collaborators: 


Participant  Type:  Postdoctoral  (scholar,  fellow  or  other  postdoctoral  position) 
Participant:  Jaedal  Jung 
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Person  Months  Worked:  4.00  Funding  Support: 

Project  Contribution: 

International  Collaboration: 

International  Travel: 

National  Academy  Member:  N 
Other  Collaborators: 


Participant  Type:  Postdoctoral  (scholar,  fellow  or  other  postdoctoral  position) 
Participant:  Mehdi  Naderi 

Person  Months  Worked:  6.00  Funding  Support: 

Project  Contribution: 

International  Collaboration: 

International  Travel: 

National  Academy  Member:  N 
Other  Collaborators: 


Participant  Type:  Graduate  Student  (research  assistant) 

Participant:  Briggitte  Morales 

Person  Months  Worked:  10.00  Funding  Support: 

Project  Contribution: 

International  Collaboration: 

International  Travel: 

National  Academy  Member:  N 
Other  Collaborators: 


Participant  Type:  Graduate  Student  (research  assistant) 

Participant:  Yuhang  Yang 

Person  Months  Worked:  12.00  Funding  Support: 

Project  Contribution: 

International  Collaboration: 

International  Travel: 

National  Academy  Member:  N 
Other  Collaborators: 


Participant  Type:  Graduate  Student  (research  assistant) 

Participant:  YunWei  Xu 

Person  Months  Worked:  12.00  Funding  Support: 

Project  Contribution: 

International  Collaboration: 

International  Travel: 

National  Academy  Member:  N 
Other  Collaborators: 


Participant  Type:  Other  (specify) 

Participant:  Scott  Floerke 

Person  Months  Worked:  10.00  Funding  Support: 

Project  Contribution: 

International  Collaboration: 

International  Travel: 

National  Academy  Member:  N 
Other  Collaborators: 


Participant  Type:  Other  (specify) 
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Participant:  James  Robertson 

Person  Months  Worked:  8.00  Funding  Support: 

Project  Contribution: 

International  Collaboration: 

International  Travel: 

National  Academy  Member:  N 
Other  Collaborators: 


ARTICLES: 

Publication  Type:  Journal  Article  Peer  Reviewed:  Y  Publication  Status:  1-Published 

Journal:  International  Journal  for  Numerical  Methods  in  Engineering 

Publication  Identifier  Type:  DOI  Publication  Identifier:  10.1002/nme.4697 

Volume:  0  Issue:  0  First  Page  #:  0 

Date  Submitted:  Date  Published: 

Publication  Location: 

Article  Title:  An  efficient  augmented  finite  element  method  for  arbitrary  cracking  and  crack  interaction  in  solids 

Authors: 

Keywords:  A-FEM,  crack  interaction,  cohesive  zone  models 

Abstract:  This  paper  presents  an  augmentation  method  that  enables  bilinear  finite  elements  to  efficiently  and 
accurately  account  for  arbitrary,  multiple  intra-elemental  discontinuities  in  2D  solids.  The  augmented  finite  element 
method  (A-FEM)  employs  four  internal  nodes  to  account  for  the  crack  displacements  due  to  an  intra-elemental 
weak  or  strong  discontinuity,  and  it  permits  repeated  elemental  augmentation  to  include  multiple  interactive 
cracks.  It  thus  enables  a  unified  treatment  of  damage  evolution  from  a  weak  discontinuity  to  a  strong  discontinuity, 
and  to  multiple  interactive  cohesive  cracks,  all  within  a  single  bilinear  element  that  employs  standard  external 
nodal  DoFs  only.  A  novel  elemental  condensation  procedure  has  been  developed  to  solve  the  internal  nodal  DoFs 
as  functions  of  the  external  nodal  DoFs  for  any  irreversible,  piece-wise  linear  cohesive  laws.  It  leads  to  a  fully 
condensed  elemental  equilibrium  equation  with  mathematical  exactness,  eliminating  the  need  for  nonlinear  equil 
Distribution  Statement:  1-Approved  for  public  release;  distribution  is  unlimited. 

Acknowledged  Federal  Support: 


Publication  Type:  Journal  Article  Peer  Reviewed:  Y  Publication  Status:  1-Published 

Journal:  Computer  Methods  in  Applied  Mechanics  and  Engineering 

Publication  Identifier  Type:  DOI  Publication  Identifier:  10.1016/j.cma.2014.03.023 

Volume:  276  Issue:  0  First  Page  #:  0 

Date  Submitted:  Date  Published: 

Publication  Location: 

Article  Title:  Generalized  bridging  domain  method  for  coupling  finite  elements  with  discrete  elements 

Authors: 

Keywords:  Discrete  element  method;  Finite  element  method;  Multiscale  method;  Spurious  wave  reflection; 
Granular  flow  phenomenon 

Abstract:  The  concurrent  coupling  of  finite  elements  and  discrete  elements  is  an  effective  DOF  reduction 
methodology  for  reproducing  the  granular  flow  phenomenon  as  discrete  element  method  does.  In  this  paper,  we 
present  a  novel  coupling  strategy  named  the  generalized  bridging  domain  method.  This  method  introduces 
independent  functions  to  weight  the  material  properties  of  the  continuum  and  those  of  the  discrete  element  model, 
and  then  the  equilibrium  of  each  new  model  is  guaranteed  by  imposing  compensation  forces.  The  compensation 
force  of  continuum  is  determined  by  force  and  displacement  compatibility  requirements  of  continuum  with  respect 
to  discrete  elements,  and  vice  versa.  Utilizing  different  weighting  functions,  four  typical  coupling  methods,  the 
bridging  domain  method,  edge-to-edge  coupling,  separate  domain  coupling,  and  separate  edge  coupling,  are 
obtained.  Additionally,  a  new  integration  algorithm  with  multiple  time  steps  is  developed  for  the  separate  edge 
coupling.  The  numerical  perfor 

Distribution  Statement:  1-Approved  for  public  release;  distribution  is  unlimited. 
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Publication  Type:  Journal  Article  Peer  Reviewed:  Y  Publication  Status:  1-Published 

Journal:  International  Journal  for  Numerical  Methods  in  Engineering 

Publication  Identifier  Type:  DOI  Publication  Identifier:  10.1002/nme.4689 

Volume:  99  Issue:  7  First  Page  #:  0 

Date  Submitted:  Date  Published: 

Publication  Location: 

Article  Title:  A  finite  element  method  with  mesh-separation-based  approximation  technique  and  its  application  in 
modeling  crack  propagation  with  adaptive  mesh  refinement 
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AUGMENTED  FINITE  ELEMENT  METHOD  FOR  HIGH-FIDELITY 
ANALYSIS  OF  STRUCTURE  COMPOSITES  (FINAL  REPORT) 

EXECUTIVE  SUMMARY 

The  proposed  objective  of  this  study,  according  to  the  original  propose,  is  to  develop  “an  efficient  and 
accurate  multiscale  computational  methodology  for  rapid  virtual  testing  of  composites  with  complex 
microstructures”.  The  proposed  research  goals  are: 

“1)  Development  of  new  2D  and  3D  A-FEM  elements  that  can  explicitly  account  for  statistical  material 
heterogeneity  and  can  faithfully  predict  progressive  damage  initiation,  propagation  and  multi-crack 
interaction”; 

“2)  Apply  the  new  method  to  study  composites  of  interest  to  the  Army,  coupled  with  3D  pCT-based 
precision  experiments,  to  understand  key  physical  phenomena  of  crack  formation  and  growth,  and 
to  quantify  their  direct  effects  on  structural  integrity  under  general  loads”; 

Two  more  goals  were  added  in  the  4*  year  extension  proposal  in  view  of  the  importance  of  material 
nonlinearity  due  to  finer  scale  (micron-  or  submicron  scale)  material  damage, 

“3)  Complete  the  development  of  a  ‘multi-physics’  nonlinear  element  within  the  AFEM  framework 
that  will  allow  us  to  represent  a  microcrack  in  a  heterogeneous  medium”; 

“4)  Verify  the  ability  of  the  multi -physics  element  to  make  correct  predictions  of  the  conditions  for 
crack  deflection,  bifurcation,  and  branching  and  of  the  influence  on  these  phenomena  of  fine-scale 
continuum  damage  in  the  surrounding  material”. 

All  these  objectives  have  been  met.  Here  we  highlight  the  key  achievements  with  more  details  given  in 
the  main  body  of  this  final  report. 

Major  Technical  Achievements  : 

•  A  new  augmentation  formulation  that  enables  standard  finite  element  method  to  deal  with  strong 
and  weak  discontinuities  caused  by  arbitrary,  coupled  multiple  damage  formation  and  propagation 
in  complex  heterogeneous,  without  the  need  for  adding  extra  DoEs  into  the  global  problem. 

•  A  novel  elemental  condensation  procedure  (consistency-check  based  algorithm)  that  demonstrated 
2-3  orders  of  magnitude  improvement  in  numerical  efficiency,  accuracy,  and  robustness. 

•  A  novel  inertia-based  stabilizing  method  that  guarantees  unconditional  convergence  and  allows  for 
smooth  transition  from  quasi-static  crack  development  to  rapid  or  even  dynamic  crack  propagation, 
and  vice  versa. 

•  A  nonlinear  A-EEM  to  account  for  both  finer  scale  (i.e.  submicron  scale)  continuum  damage  or 
plasticity  and  larger  scale  nonlinearity  due  to  nonlinear  coupling  of  multiple  cracks. 

•  Extended  the  new  A-EEM  with  coupled  thermal-mechanical  analysis  capabilities  under  static, 
cyclic,  and  dynamic  loading  conditions. 

•  Much  improved  understanding  of  the  interaction  of  multiple  damage  modes  in  composites  and  their 
convoluted  influence  on  structural  integrity  through  comparing/correlating  detailed  A-EEM  analysis 
results  against  high  resolution  experimental  measurements/characterizations  (DIC,  micro-CT,  etc) 

•  18  published/accepted  journal  publications  and  5  more  in  preparation  {Appendix  A). 

Collaborations  and  Technology  Transfer  (Appendix  C) 

•  Key  enabling  capabilities  of  the  A-EEM  have  been  enriched  and  applied  to  several  other 
governmental  and  industrial  programs. 

•  Broad  collaborations  with  nationaVinternational  universities,  research  institutions,  and  industrial 
companies. 

Graduate  Students  Involved  During  Reporting  Period  (Appendix  B) 

•  This  project  benefitted  4  PhD  students,  2  MS  students,  and  2  post-doctorate  researchers 
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1.  BACKGROUND  OVERVIEW 


Composite  materials  are  now  widely  used  in  many  key  areas  sueh  as  aerospaee,  defense, 
renewable  energy,  and  off-shore  engineering.  Reeent  high  profile  examples  include  the  Boeing 
787  dream-liners,  ultra- light  weight  composite  blades  for  wind  mills,  and  multilayer  FRPC-steel 
pipelines  for  deep-sea  oil  transportation.  Multiscale  hierarchical  microstructures  are  usually 
integrated  into  such  structures  to  achieve  multi-functional  purposes  and/or  to  further  gain  the 
weight-savings.  Key  advantages  of  composites  over  the  traditional  structural  metal  and  metal 
alloys  include  significant  weight  savings,  better  reliability  and  durability,  and  lower  manufacturing 
and  maintenance  costs.  In  particular,  composite  engineering  offers  the  exciting  possibility  to 
design  the  hierarchical  features  in  micro-  or  meso-scale  to  achieve  the  desired  structural 
performance  through  optimization  process  [1-3]. 

However,  due  to  the  relatively  short  history  of  composite  usage  as  primary  load-bearing 
structural  components,  the  current  knowledge  base  remains  inadequate  for  confident  design  of  safe 
and  durable  composite  structures.  The  safety  and  durability  uncertainty  of  composite  structures 
remains  high  and  damage  tolerance  analysis  must  be  considered  as  an  integral  part  in  the  composite 
structural  design  process.  However,  unlike  structural  metals  which  are  homogeneous  and  isotropic, 
composites  are  inherently  inhomogeneous  and  anisotropic.  Microscopic  flaws  and  imperfections 
are  inevitably  present  in  these  materials  due  to  manufacturing  processes.  Currently  the  damage 
tolerance  and  durability  design  relies  heavily  on  large  experimental  test  programs  that  are 
extremely  costly  and  time  consuming.  High-fidelity  simulation  capabilities  that  can  potentially 
offer  efficient  and  accurate  strength  and  durability  predictions  based  on  constituent  materials 
properties  can  greatly  reduce  the  costly  experimental  inputs.  However,  accurate 
estimation/prediction  of  composite  strength  at  structural  scale  remains  elusive  despite  decades  of 
intensive  research  [4-6]. 

It  proves  very  changing  to  faithfully  predict  the  ultimate  strength  of  composites  due  to  the 
fact  that  almost  all  composites  exhibit  complex  progressive  damage  evolution  phenomena  before 
complete  failure.  For  example,  a  composite  laminate  will  typically  develop  non-critical  intra-ply 
cracks  (transverse  cracks  in  off-axial  plies  and  longitudinal  splits  in  0°-plies)  far  below  the  ultimate 
load,  which  facilitate  local  delaminations  between  neighboring  plies.  These  damage  processes 
continuously  evolve  in  a  nonlinearly  coupled  manner  and  eventually  lead  to  the  formation  of  large 
delaminations,  which  in  turn  cause  significant  fiber  breakage  in  primary  load  bearing  plies  and 
lead  to  ultimate  stmctural  failure  [7,  8].  While  it  is  desired  from  damage  tolerance  point  of  view, 
predicting  such  progressive  damage  evolution  up  to  ultimate  failure  proves  to  be  very  challenging. 

Thanks  to  recent  development  in  numerical  modeling  methods  and  high-resolution 
experimental  characterization,  significant  advances  in  composite  strength  prediction,  especially  in 
continuous  fiber  reinforced  composite  laminates,  have  been  made  in  the  past  decade  [1,  2,  9,  10]. 
This  was  enabled  by  1)  the  development  of  advanced  finite  element  methods  (FEMs)  including 
the  extended  finite  element  method  (X-FEM,  e.g.  [11-14]),  Phantom  node  method  (PNM)  (e.g. 
[15-17]),  and  the  PNM -based  augmented  finite  element  method  (PNM-based  A-EEM,  e.g.  [18-21], 
and  2)  the  nonlinear  cohesive  zone  models  (CZMs)  for  discrete  fracture  events  in  composites  (e.g. 
[22-28].  The  key  is  to  embed  CZM  descriptions  for  major  damage  modes  into  the  advanced  EEM 
formulations  so  that  the  discrete  cracking  events  can  be  initiated  and  propagated  at  the  locations 
determined  by  local  stresses  in  response  to  external  loads.  Eor  a  more  comprehensive  review  of 
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the  recent  advances  and  the  remaining  outstanding  numerical  and  material  issues  associated  with 
such  progressive  damage  modeling  in  composite  laminates,  the  readers  are  referred  to  the  recent 
review  article  by  Van  de  Meer  [9]. 

It  has  been  generally  recognized  that,  to  achieve  high-fidelity  predictions  of  the  progressive 
failure  processes  in  structural  composites,  several  key  elements  have  to  be  considered 
simultaneously  in  the  numerical  framework:  1)  a  numerical  model  should  allow  for 
crack/discontinuity  initiation  from  a  pristine  continuum  at  locations  determined  by  local  stress 
environments;  2)  carefully  calibrated  cohesive  laws  for  crack/discontinuity  evolution  and 
propagation;  3)  crack/discontinuity  initiation  criteria  that  account  for  the  constraining  effects  from 
neighboring  plies  (in-situ  strength,  e.g.  [29,  30]);  and  4)  correct  representation  of  the  shear 
nonlinearity  at  ply  level  [20,  25,  3 1];  4)  special  numerical  treatment  is  needed  to  correctly  account 
for  the  coupling  between  the  cracked  ply  elements  and  the  interfacial  cohesive  elements  (otherwise 
correct  stress  transfer  is  not  always  guaranteed  and  numerical  stress  locking  may  occur  [32-34]); 
and  5)  a  computationally  efficient  algorithm  that  can  track  all  individual  cracks  and  their  nonlinear 
coupling,  preferably  without  the  need  to  introduce  additional  DoFs  (The  additional  DoFs  in  X- 
FEM  and  PNM-FEM  are  specific  to  a  single  crack  and  they  become  problematic  when  multiple 
crack  merge  or  a  single  crack  branches  into  multiple  cracks). 

Eurther,  in  textile  composites  the  material  heterogeneity  issue  of  composites  cannot  be 
resolved  adequately  by  mainstream  formulations  of  conventional  materials/structures  modeling, 
owing  to  the  complex  interaction  between  fiber  reinforcement  architecture  and  constituent  material 
heterogeneity,  which  are  on  the  same  scale  as  that  of  the  features  of  the  structures.  This  negates 
the  common  strategy  of  homogenizing  the  material  properties  in  simulations.  The  heterogeneity 
arises  from  the  distinct  fiber  tows  in  the  textile,  which  are  strongly  anisotropic  and  their  shapes 
may  be  significantly  distorted  due  to  manufacturing  process.  The  tows  are  typically  ~  0.1  -  1  mm 
in  lateral  dimension,  which  is  of  the  same  order  of  the  panel  thickness  (1~  10  mm).  Eurthermore, 
the  interlacing  pattern  of  an  integral  textile  structure,  which  has  profound  impact  on  the  elastic 
properties  and  failure  strength  [35-37],  can  be  quite  complex,  involving  many  hundreds  of 
geometrically  dissimilar  tow  segments  or  local  tow  arrangements.  The  structure  will  generally  be 
subject  to  non-periodic  mechanical  and  thermal  loads,  including  spatial  and  temporal  peaks.  In  a 
piece  of  structure  of  relevant  size,  it  remains  a  serious  challenge  to  develop  a  mesh  that  deals  with 
the  interlacing  architecture,  the  heterogeneity  of  the  tows,  and  general  loads,  while  maintaining 
computational  tractability.  Strength  prediction  is  even  more  challenging:  the  randomly  distributed 
matrix-rich  pockets  have  to  be  explicitly  considered  because  they  are  susceptible  to  microcrack 
formation  at  relatively  small  loads.  The  complex  coupling  of  all  the  nonlinear  processes  associated 
with  crack  nucleation,  bifurcation  and  deflection  must  be  adequately  captured. 

The  above  challenges  require  the  establishment  of  numerical  modeling  platforms  that  can 
deal  with  not  only  the  arbitrary  cracking  issue,  but  also  the  issue  of  arbitrary  crack  interaction  such 
as  crack  merging  and  crack  bifurcation,  again,  without  any  a  priori  knowledge  of  the  locations  and 
load  levels  at  which  such  events  should  happen.  Unfortunately,  the  above  mentioned  advanced 
numerical  methods  (X-EEM,  PNM-based  A-EEM)  are  all  very  inflexible  in  dealing  with  multiple 
crack  merging  or  bifurcation.  This  is  because,  in  these  methods  each  individual  crack  requires  its 
own  additional  copy  of  degree  of  freedoms  (DoEs)  (X-EEM)  or  extra  node  copies  (PNM-based  A- 
EEM)  to  describe  the  discontinuous  displacement  field  associated  with  it.  Tracing  the  evolution  of 
complex  fracture  surfaces  of  multiple  cracks  quickly  becomes  extremely  tedious  and  numerically 
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burdensome  because,  in  such  cases,  there  are  no  effective  tracking  algorithms  that  can  judiciously 
assign  different  copies  of  DoFs  or  nodes  to  different  cracks  for  complex  crack  configurations. 

It  is  therefore  beneficial  to  seek  numerical  methods  that  can  account  for  arbitrary 
discontinuities  and  their  nonlinear  interactions  with  less  numerical  burden  for  tracing  complex 
crack  surfaces.  In  this  regard,  there  are  mainly  two  methods  to  reduce  or  eliminate  the  need  for 
additional  DoFs  or  nodes.  One  method  is  the  phase-field  method  [38-40].  This  method  introduces 
an  additional  nodal  DoF  to  approximate  a  fracture  surface  with  a  highly  concentrated  yet 
continuous  phase-field,  which  smoothes  the  boundary  of  the  crack  over  a  small  region.  The  major 
advantage  is  that  the  evolution  of  the  fracture  surfaces  follows  from  the  solution  of  a  coupled 
system  of  partial  differential  equations.  Thus,  it  does  not  require  the  fracture  surfaces  to  be  tracked 
algorithmically.  However,  this  method  is  mesh-dependent  and  it  requires  an  extremely  fine  mesh 
to  resolve  the  sharp  discontinuity  associated  with  an  interface  or  a  crack  surface.  The 
computational  cost  remains  extremely  high.  There  are  also  serious  ambiguities  regarding  the  non¬ 
zero  stresses  at  crack  wake  surfaces. 

Another  approach  is  the  embedded  discontinuity  method.  This  method  follows  the  earlier 
work  in  smeared  localization  models  and  enhanced  strain  methods  [41-46]  and  it  has  been 
extended  to  account  for  discrete  discontinuities  by  several  research  groups  [47-54]  following  Simo 
and  colleagues’  seminal  approach  in  treating  plastic  strain  localizations  with  discrete  discontinuous 
descriptions  [55,  56].  Special  (or  regularized)  shape  functions  are  used  to  account  for  the 
discontinuous  crack  displacements  within  an  element  [50,  51,  53,  56].  However,  the  special  shape 
functions  have  to  be  carefully  constructed  so  that  the  orthogonal  principle  is  satisfied.  Otherwise, 
spurious  deformation  modes  or  serious  stress  locking  may  occur  [54,  55].  Furthermore,  the  ability 
of  this  approach  in  handling  interactive  cracks  is  limited  —  thus  far  there  is  no  report  on  regularized 
shape  functions  that  can  descript  multiple  discontinuities  within  a  single  element. 

Thus  the  main  focus  of  this  study  is  to  develop  new  finite  element  methods  that  can  deal 
with  discrete  cracks  without  the  need  of  introducing  additional  DoFs  or  nodes,  so  that  the 
numerical  burden  for  tracing  crack  surfaces  can  be  minimized.  During  the  proposed  research 
period,  we  have  formulated  and  implemented  a  novel  A-FEM  formulation  that  can  treat  multiple 
intra-element  cohesive  cracks  in  homogeneous  or  heterogeneous  solids,  without  the  need  to 
introduce  any  extra  DoFs.  In  the  new  A-FEM  formulation  four  internal  nodes  with  regular 
displacement  DoEs  were  introduced  to  account  for  a  cohesive  crack,  so  that  the  crack 
displacements  are  natural  outcomes  from  element  equilibrium  consideration.  Thus  it  does  not  need 
to  assume  the  deformation  modes  a  priori  for  crack  displacements  as  in  the  enhanced  strain  method 
or  in  many  of  the  embedded  discontinuity  methods.  More  importantly,  we  demonstrated  the  new 
A-EEM  can  account  for  1)  intra-element  material  heterogeneity  which  involves  a  natural  evolution 
from  a  weak  discontinuity  to  a  strong  one,  and  2)  repeated  elemental  augmentation  to  enable 
multiple,  interactive  intra-element  discontinuities  to  account  for  crack  merging  (from  a  bulk 
material  domain  into  a  material  interface)  or  crack  bifurcation  away  from  the  interface.  In  addition, 
we  report  a  novel  and  efficient  solving  algorithm  for  elemental  condensation  that  provides 
analytical  solutions  to  the  local  equilibrium  equations  with  embedded  piece- wise  linear  cohesive 
crack  like  discontinuities.  It  will  be  demonstrated  that  the  new  A-EEM  with  the  help  of  the  new 
condensation  algorithm  can  provide  orders  of  magnitude  improvement  in  numerical  efficiency, 
accuracy  and  robustness  [21,  57]. 
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2.  PROBLEM  STATEMENT 


2.1  Strong  and  Weak  Form 

Consider  a  discontinuous  physical  domain  Q  as  shown  in  Figure  1.  The  domain  is 
composed  of  two  sub-domains  of  the  same  or  different  materials,  i.  e.,  Q  =  uQ  .  The  two 
domains  are  connected  by  a  discontinuity  which  could  be  a  cohesive  crack  or  a  material  boundary, 
r.  =  Cur/.>  The  prescribed  tractions  /  are  imposed  on  boundary  Ff  and  the  prescribed 

displacements  m  on  boundary  Fu.  The  stress  field  inside  the  domain,  ct,  is  related  to  the  external 
loading/.  The  strong  form  of  the  equilibrium  equations  and  boundary  conditions  are: 


V-(T^=0  (inQ^); 

rs  n  =r  (onF/) 
u^=r  (onF/) 
t  =  a  n  =-t  (onF/) 


V  •  o  =0  (in  Q  ) 

(T  -n"  =  f'  (on  F^") 
u’=u'  (onF/) 
t  =(T~-n"  =  t  (onF/) 


(1) 


where  and  a  are  the  stresses  in  the  subdomains  and  Q  ,  respectively;  and  t 

are  the  tractions  along  the  discontinuity  surfaces  F^"^  and  F^” ,  respectively;  and  are  the 

outward  normal  of  F^^  and  F^" ,  respectively.  The  last  two  expressions  in  Eqn  (1)  come  from  the 
stress  continuity  across  the  discontinuity.  The  traction,  t,  is  a  function  of  the  relative  displacements, 
Au  ,  between  F^^  and  F^” ,  i.e.. 


t  =  t(Au)  (onFJ 

Au  =  u^-u^  (onFJ 


where  u'^  and  u'  are  the  displacement  fields  in  and  respectively.  Eqn  (2)  serves  as  the 
constitutive  law  and  kinematic  equation  of  the  discontinuity.  In  this  study,  cohesive  zone  models 
with  piece-wise  linear  traction- separation  laws  were  used  to  describe  strong  discontinuities. 


Figure  1 :  Notation  for  an  elastic  body  with  a  discontinuity 


'  For  a  weak  discontinuity  while  for  a  strong  discontinuity  and  are  separated. 


5 


The  domain  surrounding  the  discontinuity  is  assumed  to  be  elastic.  We  further  assume 
small  strains  and  displacements  condition.  Thus  the  constitutive  law  and  kinematic  equations  for 
the  two  domains  can  be  written  as: 

(inQ^);  (inQ’) 

£+  -  [Vu^  +  (Vu^)'^]/2  (in  QO;  s  -  [Vu'  +  (Vu')^]/2  (in  Q^) 

where  and  C”  are  the  material  stiffness  tensors  of  the  two  sub-domains  traversed  by  the 
discontinuity,  respectively.  They  can  be  of  isotropic  or  orthotropic.  The  superscript  (O^  denotes 
transposition. 

The  above  strong  form  of  the  problem  can  be  converted  into  a  weak  form  using  the 
principle  of  virtual  work.  The  displacement  fields  (  and  u  )  must  be  sub-sets  of  the 
kinematically  admissible  displacement  fields,  IJ: 

u""  elJ  =  {v^  e  V;v^  =u^  onT^""};  u“  elJ  =  {v“  eV;v”  =u“  onT^”}.  (4) 

Applying  the  principle  of  virtual  work  separately  to  the  subdomains,  the  weak  forms  can  be  written 
as 


[  :  £^(v^)  dQ  =  f  F^-v^dT-f  t(Av)-v^dr 

J  Cl  J  c  p-  J  r 

f  tr  :£(v)dQ=f  F  v  dT-i-f  t(Av)  v  dT 

J  J  Fr-"  J 


V(v^V■)GlJ 


(5) 


where  Av  =  v'^-v~  (onT^) .  The  left-hand-sides  are  the  internal  virtual  work  and  the  right- 
hand-sides  are  the  virtual  work  done  by  the  external  forces  and  the  tractions  along  the  crack 
surfaces.  Note  that,  if  one  of  the  subdomains  (say  QT )  further  cracks,  the  principle  of  virtual  work 
can  be  applied  to  all  three  subdomains  without  any  further  modifications. 

2.2.  Piece-wise  Linear  Cohesive  Laws 

The  discontinuity  in  section  2.1  is  often  modeled  by  cohesive  laws  that  provide  explicit 
traction-separation  relations.  Any  nonlinear  cohesive  laws  can  be  approximated  with  piece-wise 
linear  laws  as  shown  in  Figure  2.  The  cohesive  law  employs  separate  mode  I  and  mode  II  traction- 
separation  relations  and  it  was  first  proposed  by  Yang,  Thouless  and  colleagues  [22,  58-60].  For 
those  cohesive  laws  employing  coupled  mode  I  and  mode  II  relation  as  in  [24,  26],  the  piece-wise 
linear  form  of  the  traction-separation  laws  can  be  obtained  with  an  incremental  form. 
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Figure  2  A  mixed  mode  Piece-wise  linear  cohesive  law  with  (a)  mode  I  and  (b)  mode  II  traction  separation 
relations. 

For  the  piece- wise  linearized  models,  it  is  advantageous  to  index  the  traction-separation 
laws  by  cohesive  segment  numbers  as  labeled  in  Figure  2  by  the  circled  numbers.  Each  segment 
is  characterized  by  a  critical  stress  (f^'^  or  ),  a  corresponding  local  crack  displacement 

or  5^^'' ),  and  a  constant  stiffness  ( or  ).  Here  the  superscripts  (/)  and  (j)  (i,  j  =  1,2,  and 

3)  are  free  indices  for  the  cohesive  segments  with  constant  stiffness.  For  the  cohesive  model  shown 
in  Figure  2, 


(T**’  =  o^^;  =  (Tj;  = 


0; 


P  *'  ‘•2’ 

=  A  •  ^  e  .  e(3)  e 

‘^nP  ^ sV  ^ sl^  ^ sc  ■ 


(6) 

(7) 


And  the  constant  cohesive  slopes  for  any  segments  are: 

c,;')  =  (f« {i,j  =  1,2,3) ,  (8) 

where  6-^°^  =  =  0  and  -  0 .  Further,  the  crack  displacement  ranges  consistent  with 

the  cohesive  segments  are  indexed  with 


(b;  =  i,2,3) 


(9) 


Thus  the  linearized  relation  between  the  cohesive  stresses,  and  t(SJ,  and  the  crack 

displacements,  5^  and  5^,  for  any  segments  may  be  written  as: 


t{5^  )  =  sgn((^^ )  -  (J/'-'’ )]  (for  e  5s'^  &&  |  >  ) 


(7-1)' 


(;■) 


(10) 


(for  S^eSn  &&  ) 
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where  the  sign  function  sgn(»)  is  defined  as  sgn(»)  =  < 


if  •  >  0^ 
if  *  =  0 
if  *<0 


denotes  taking  the 


1 

0 

-1 


absolute  value  of  the  argument.  and  are  two  solution-dependent  variables  related  to  the 
irreversibility  of  the  cohesive  model,  as  discussed  below. 


The  irreversibility  of  the  cohesive  model  was  explicitly  considered  by  an  additional 
unloading  segment  in  each  traction- separation  law  (segment  4  in  Figure  2).  In  this  paper,  the 
maximum  crack  displacements  ever  reached,  <5'„*  for  mode  I  and  <5’^,  for  mode  II  as  shown  in 
Figure  2,  were  used  as  the  historical  (or,  solution-dependent)  variables  to  facilitate  distinguishing 
between  loading  and  unloading  path,  i.e.,  unloading  occurs  when  the  current  crack  displacements 
are  smaller  than  the  respective  historical  displacements.  The  corresponding  historical  cohesive 
stresses,  <7*  and  f * ,  were  computed  from  these  historical  crack  displacements  using  Eqn  (10). 
Thus  for  segment  i  or  j  =  4,  the  cohesive  stress  -  crack  displacement  relation  is 


(4) 

)  =  sgn((J^  )«f  ’  ;  (for  \S^\eS,  &&  ) 

ct(S^  )  =  (for  e  S„  && 


The  cohesive  stiffnesses  and  the  displacement  ranges  are 


(12) 

v;  (4) 

v;  (4) 

Ss  =[0,<^^*]; 

dn  =[-CO,<5^^J 

For  other  cohesive  laws  such  cohesive  segments  may  be  constructed  similarly  after  piece- 
wise  linearization. 

Based  on  the  individual  modes,  a  mixed-mode  cohesive  model  can  be  constructed  by 
recognizing  that  the  total  traction- separation  work  absorbed  during  fracture,  can  be  separated 
into  the  mode  I  and  mode  II  components,  and  ,  so  that. 


(13) 


The  two  separate  components  can  be  calculated  by  integration  of  the  mode-I  and  mode-II  traction- 
separation  curves  (Fig.  2): 


g,=J(V(<5)d5 

5„=|(‘r(5)d5 


(14) 


Note  that  Sn  and  &  are  not  independent  parameters;  they  evolve  together  as  a  natural  result 
of  the  interplay  between  the  deformation  of  two  joined  domains  and  the  details  of  the  two  pure- 


8 


mode  traction- separation  laws.  A  failure  criterion  is  required  to  determine  the  critical  values  of 
the  two  components  of  (^,  and  (shaded  areas  in  Figure  2),  at  which  complete  fracture  of 
the  cohesive  zone  occurs.  The  criterion  used  in  this  study  is  a  simple  one  from  [61] 

g;/r,+g;/r„  =  i  (i5) 

where  Fi  and  Fii  are  the  total  areas  under  the  pure  mode  I  and  pure  mode  II  cohesive  laws.  They 
are  the  mode  I  and  mode  II  fracture  toughnesses  in  the  linear-elastic-fracture-mechanics  (LEFM) 
context. 

A  more  detailed  account  of  this  mixed-mode  cohesive  zone  model  can  be  found  in  [23,  59, 
60].  The  major  advantage  of  this  cohesive  law  is  that  there  is  no  need  to  specify  the  mode 
mixedness  a  priori.  The  mode  mixedness  and  the  mixed-mode  toughness  evolve  as  numerical 
outcomes  of  the  local  equilibrium  of  stresses.  More  importantly,  this  law  guarantees  correct  mode 
mixedness  when  LEFM  conditions  are  satisfied  [34,  62,  63]. 


3.  A  NEW  2-D  A-FEM  FORMULATION 


3.1  A-FEM  Formulation  with  Single  Intra-Elemental  Crack 

Without  the  loss  of  generality,  the  4-node,  quadrilateral  plane  element  is  used  here  to 
illustrate  the  augmented  finite  element  scheme.  The  physical  element  with  regular  or  external 
nodes  1,  2,  3,  and  4  is  severed  by  a  cohesive  crack  or  a  bi-material  interface  as  shown  in  Figure  3. 
We  take  a  common  approach  that  assumes:  1)  during  crack  growth,  a  crack  tip  always  resides  at 
an  element  boundary  during  its  propagation  as  in  the  phantom  node  method  and  the  alike  [18,  20, 
64,  65],  and  2)  the  existence  of  a  cohesive  zone  eliminates  the  crack-tip  singularity  [22,  23,  66- 
68].  Thus  the  method  proposed  here  cannot  treat  cracks  with  crack-tip  singularities  as  can  be 
conveniently  done  in  the  X-FEM  framework.  However,  for  many  engineering  materials  that  fail 
by  the  development  of  finite  damage  process  zones  (i.e.,  polymers,  concretes,  ductile  metals,  fiber- 
reinforced  polymeric  composites,  etc),  such  a  treatment  is  appropriate  [23,  28,  65,  69,  70]. 


If  the  element  is  cut  by  a  single  crack,  there  are  two  basic  crack  configurations:  1)  two 
quadrilateral  subdomains  as  shown  in  Figure  3(b),  and  2)  one  triangular  subdomain  and  one 
pentagonal  subdomain  as  shown  in  Figure  (3c).  For  each  crack,  a  local  coordinate  system  is 
defined  by  the  directions  along  (s)  and  perpendicular  to  (n)  the  crack  path,  from  which  the  normal 
(5n)  and  shear  (5s)  crack  displacements  can  be  defined  as  shown  in  Figure  3(b).  These  local  crack 
displacements  are  related  to  the  global  crack  displacements.  Am  and  Av,  by 


cos  sin 


is  the  rotational  matrix  between  the 


where  [RJ-  . 

[_-sm6'i  cos  O'! 

global  and  local  coordinates  and  6\  is  the  angle  between  global  x-axis  and  local  s-axis  (Figure  3). 


We  introduce  4  internal  nodes  (node  5,  6,  6’,  and  5’  in  Figure  3  (b)  and  (c))  with  regular 
displacement  DoFs  so  that  the  displacement  jumps  across  the  crack  are  simply  the  differences  in 
displacements  between  the  respective  node -pairs  5-5’  and  6-6’. 
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With  the  help  of  these  internal  nodes,  the  displaeement  field  in  eaeh  subdomain, 
u“  (a  =  +  or  -) ,  ean  be  obtained  through  standard  FE  shape  funetion  interpolation  as  below, 

u“=N“-d“;  (16) 


where  d"  is  the  nodal  displacement  array  of  the  subdomain  a  (cr  =  +  or  -) .  If  the  element  is  cut 
into  2  quadrilateral  sub-domains  as  shown  in  Figure  3(b),  d^  =  {u^,v^,u^,v^,\u^,,v^.,Uy,v^Y  and 

d~  =  {Mj,Vj,M2,V2,I  Mg,Vg,M5,V5}^ .  N“  is  the  interpolation  matrix  composed  of  standard  bi-linear 
isoparametric  shape  functions.  The  symbol  “|”  indicates  that  the  displacements  ahead  of  it  are 
external  nodal  displacements  which  are  denoted  with  ,  and  those  behind  it  are 

the  internal  nodal  displacements  which  are  denoted  with  .  Similarly, 

d;,  ={Mi,Vi,a2,V2}^  and  . 


Figure  3  Illustration  of  the  element  augmentation  from  (a)  a  regular  element  with  possible  different  material 
domains,  to  (b)  an  A-FE  with  two  quadrilateral  sub-domains,  or  to  (c)  an  A-FE  with  one  triangular  sub-domain 
and  one  pentagonal  sub-domain. 


If  the  element  is  cut  into  a  triangular  and  a  pentagonal  subdomain  as  in  Figure  3(c),  then 


fd\l 

1 

\d-,] 

J  ext  \ 

and  d  =j 

1  ext  [ 

1  ,  f 

iCj 

KtJ 

{  Mp  V| ,  1^2  ?  F2 , 1^2 ,  ,  I  } 


The 


shape  functions  in  the  interpolation  matrix  N“  are  from  the  polygon  FEM  with  optimized 
quadrature  rules  recently  developed  by  Sukumar  and  colleagues  [71,  72].  The  strain  field  in  each 
subdomain  is: 


£“=(a^''N“)-d“=B“-d“; 


(17) 


where  is  the  differentiation  operator,  i.  e.,  d^  = 


d/dx 

0 


called  strain  matrix. 


0  d/dy 
d/dy  d/dx 


and  B“  is  the  so¬ 


lo 


Substituting  Eqn  (16)  and  (17)  into  Eqn  (5)  and  using  the  principle  of  virtual  work,  one 
obtains  the  weak  form  of  the  equilibrium  equation  for  standard  elements 


jr-d" 

[L'-d” 


=  F^  +F"" 

ext  int 


F  +F 

^  ext  '  ^  int 


(18) 


where 


If  =  r  ^(B“)'^D“B“dQ  is  the  stiffness  matrix  corresponding  to  subdomain  a; 

J  O 


F“j  =  [  „  (N“)^f  “dF  is  the  equivalent  external  force  array  on  subdomain  a; 


Fj“  =  I  „  dF  is  the  equivalent  internal  force  array  integrated  from  the  cohesive 

stresses  along  the  intra-element  cohesive  crack. 

D“  is  the  material  stiffness  matrix  transformed  from  the  local  coordinates  into  the  global 
coordinates  (see,  for  example  [73]). 

In  the  interest  of  developing  an  efficient  A-FE  that  can  deal  with  multiple  interactive  intra- 
elemental  cracks,  here  we  adopt  the  non-conforming  approach,  which  basically  assumes  that  the 
inter-element  continuity  is  enforced  only  through  the  regular  (or  external)  nodes  and  their  nodal 

DoFs  (d^^j  and  d~^,  );  while  the  conformity  of  the  internal  DoFs  associated  with  the  internal  nodes 
(dj^„,  and  dj'j,,  )  are  not  enforced.  This  simplification  leads  to  the  conclusion  that  F^"^  is  a  function 
of  regular  nodal  DoFs  and  d“^,  )  only,  and  F7  is  a  function  of  internal  DoFs  (dj^„j  and  d,”„,  ) 
only,  which  allow  us  to  re-write  Eqn  (8)  as  follows. 


(19) 


Where  L“  (i,  j  =  1,  2;  a  =  -i-,  -)  are  the  sub-matrices  of  L“ ,  i.e.,  L“  = 


T  a  T  a 

^11  ^12 

T  a  T  a 

.^21  ^22. 


.  We  note  that 


L22  is  the  sub-stiffness  matrix  related  to  the  two  internal  nodes  of  the  subdomain  a. 


Next  we  derive  the  condensed  elemental  equilibrium  equation  for  both  of  the 
configurations  in  Fig  3(b)  and  Fig  3(c)  in  a  unified  formulation.  It  is  noted  that  for  both 
configurations,  the  internal  DoF  and  equivalent  force  arrays  are  the  same,  i.  e.. 
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(20a) 


*^int  *^6'5'  }  ’ 

*^int  “  *^65  “  {*^6’^6’^5’^5  }  ’ 


F^=F  =[F  F  F  F 

^int  ^6'5'  ;t6'’^  y6'’^  ;t5'’^  y5' J 

F^=F=fF  F  F  F 

^int  ^65  ;t6’^  y6’^  ;t5’^  y5  J 


The  external  DoF  and  force  arrays  are  specific  to  a  particular  crack  configuration.  For  the  two 
quadrilateral  subdomain  configurations  (Fig  3b), 


C.r  =d34  ={M3’^3’M4’^4r; 
d;,,  =di2  ={Mi,Vi,M2>V2r; 


^ext  ~  ^34  “  { -^x3  ’  -^y3  ’  -^x4  ’  -^y4  I 

F.;,  =  ^n={Fxl^FyX^Fx2^Fy2V  ’ 


(20b) 


For  the  triangular-pentagonal  subdomain  configuration  of  Fig  3(c), 


dgxr  d4  , 

d^xr  =di23  ={Mi,Vi,M2,V2,M3,V3}'^; 


F^  =F  =fF  F  1^ 

^ext  ^4  I4x4’^.v4i 

F  =F  =fF  F  F  F  F  F 

^ext  ^123  t^  xl’^  >1’^  x2’^  >’2’^  x3’^  y3  J 


(20c) 


However,  for  both  configurations,  the  equilibrium  can  be  written  as 


Lij] 

_F2i] 


(a) 

(b) 

(c) 

(d) 


(for  Q^) 
(for  Q^) 


(21) 


Also,  due  to  stress  continuity  across  the  cohesive  crack,  Fg5=-Fg,g,  =  dF  .  Eqn  (21)  is 

valid  for  both  weak  and  strong  discontinuous  problems. 

Weak  Discontinuity:  For  weak  discontinuity  problems,  dg.5,  =  dgg .  Together  with  the  condition 
of  Fg5=  -  Fg.g, ,  Eqn  (21)  further  simplifies  to 


Fii  Lj2A  Lji  Fj2A  Lji  jdgj-J 

^  -L)2A"t2,  L;-L;2A“‘r2jld:J 


(22) 


where  A  =  L22  +  L22 .  Eurthermore,  the  internal  forces  are  determined  as: 


Fgg  ={(I-L22A-‘)l-„  -L22A-*L;}pj  (23) 

Thus,  for  an  element  containing  two  material  domains  with  a  perfectly  bonded  interface, 
once  the  external  nodal  displacements  are  given  (which  is  the  case  for  displacement-based  EEs), 
the  resultant  forces  acting  along  the  interface  are  immediately  known.  This  solution  is 
mathematically  rigorous  in  the  EEM  sense. 
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With  Fgj  known,  the  current  tractions  acting  on  the  interface  are  known.  For  example,  the 
normal  and  shear  stresses  averaged  along  the  interface,  in  local  coordinates,  are 


-  =[RJ 


,  where  is  the  interface  length  (Fig  3a).  These  stresses  can  then 


be  used  to  determine  the  initiation  of  interface  debonding,  which  is  the  turning  point  from  a  weak 
discontinuity  to  a  strong  one.  A  common  initiation  criterion  used  in  this  paper,  is 


(24) 


where  the  MacAuley  bracket  <  •  >  ,  defined  as  <•>=(•  + 1*1)  /  2 ,  is  used  to  reflect  the  physical 
observation  that  compressive  normal  stresses  do  not  contribute  to  crack  initiation. 

If  this  criterion  of  Eqn  (24)  is  met,  the  discontinuity  becomes  a  strong  one,  i.e.,  dg.j,  ^  dgg , 
which  is  addressed  next. 


Strong  Discontinuity:  For  strong  discontinuities,  we  shall  consider  cohesive  cracks  with  mixed¬ 
mode  piece- wise  linear  traction- separation  laws  as  detailed  in  Figure  2,  where  the  linear  segments 
of  each  traction- separation  laws  are  properly  indexed.  The  indexing  scheme  is  closely  related  to 
the  new  condensation  solver  that  will  be  described  shortly. 


Following  the  notation  established  in  section  2.2,  the  cohesive  stresses  between  the  two  pairs  of 
internal  nodes  5-5’  and  6-6’  can  be  expressed  in  its  local  coordinates  as 


(0 


U)  J. 

6'6 


(k) 
5'5 


f  =  o.  +  a.  }"  (i,  i, t,  /  =  1, 2, . . 4) 


(25a) 


where 


;(M) 


(7 


(/-I) 


-a. 


(25b) 


is  the  characteristic  stress  matrix.  The  critical  stresses,  and  (7^^^  (/,  7  =  1,..,4)  have  been 
defined  in  section  2.2  (Eqn  10,  Eng  1 1  and  Eigure  2). 

tto  =Diag[«/'^ ;  ]  (25c) 


is  the  cohesive  stiffness  matrix.  Here  Diag[a;  b',  c;  d]  means  a  diagonal  matrix  with  non-zero 
arguments  a,  b,  c,  and  d  at  the  corresponding  diagonal  places.  This  notation  will  be  used  in  the  rest 
of  this  paper.  is  crack  displacement  array  at  internal  node -pairs  6’-6 

and  5 ’-5  measured  in  local  coordinate  system. 

Each  of  the  free  indices  (/,  j,  k,  1)  can  range  from  1  to  the  maximum  cohesive  segment 
number,  indicating  in  which  segments  of  the  cohesive  laws  the  current  cohesive  stresses  are  located 
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(Section  2.2).  Note  that  the  characteristic  stress  (Oj)  and  stiffness  matrices  (a^)  are  determined  if 
any  particular  combination  of  (i,j,  k,  1)  is  chosen. 

Fg5  and  Fg.j,  are  obtained  by  integrating  the  cohesive  stresses  in  Eqn  (25a)  over  the  crack 
plane  and  then  converted  into  global  coordinates  using  the  rotation  matrix  Ri ,  which  gives: 

Ffis  =  -K5'  =  L  (So  +  a  •  )  (26a) 

where  -  dg.j,  -dgj  -  {(dg.  -dg),  (dj.  -dj)}^  are  the  crack  displacements  at  node  pairs  6 ’-6  and 
5 ’-5  in  global  coordinates,  i.  e., 

Ad  in,  =  dg.g,  -  dgj  =  Diag[R, ;  R,  ]  { }  (26b) 


{Sq}  and  [a]  are  related  to  the  characteristic  cohesive  stress  and  stiffness  matrices  in  Eqn  (25): 


Sg=Diag[R,^R,^](T,nJ^  a, 
a  =  Diag[R,'';  R,'"]  (T^^g)^  a,  Diag[R,;  R,] 


Here  T^nh  ^coh  the  integration  and  interpolation  matrices  associated  with  the  cohesive  stress 
integration  and  the  displacement  interpolation  along  the  crack  plane,  respectively.  Eor  the  single 
crack  case,  we  used  the  standard  Gaussian  interpolation  and  integration  scheme  with  2  integration 
points,  i.e.. 


N 


coh 


■Diag[4;4] 

Diag[4;4r 

•  T 

_k 

■Diag[4;4] 

Diag[4;4]‘ 

_Diag[4;4] 

Diag[4;4]_ 

’  ■“‘coh 

2 

_Diag[4;4] 

Diag[4;4]_ 

where  =(^l-l/V3j/2  and  =(^l  +  l/^/3j/2  .  (The  matrices  for  several  other  popular 
interpolation  and  integration  schemes  have  been  summarized  and  evaluated  in  [57,  74]). 

Substituting  Eqn  (26)  into  expressions  (b)  and  (d)  in  Eqn  (21) ,  the  following  equation  is  derived: 

J (L22  +  4i®)^6'5'  ~  ~  ~  (21) 

1~(4i®)^6'5'  (^22  4i®)^65  ~  hl^O  ~ 


Eor  a  given  cohesive  segment  combination  ((i,J,  k,  1)  this  is  a  linear  equation,  from  which  the  crack 
displacements  at  the  internal  node  pairs,  Adj^,  =  dg.j.  -dgj ,  is  solved  explicitly: 


Ad„,  A-'  +B-'  +/„B-'a(<F;)-‘)S. 

+(B-‘-;„A-y>F;,r')L;,d„,+(/„B-'o('p;,)-‘-A-‘)L;,d;„ 


(28a) 
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where, 


=  K2  +  4i«;'P22  =  L22  +  4i«;  A  =  n  -  B  =  (28b) 

The  explicit  solution  to  the  crack  displacements  can  also  be  solved  as: 


0 


0 


-B  L,, 


-a-'l;, 


-4,B-'a(T;,)-^L- 


-1t  + 
21 


(28c) 


The  fully  condensed  equilibrium  equation  without  any  internal  forces  and  displacements 
can  be  obtained  by  substituting  Eqn  (28)  into  expressions  (a)  and  (c)  in  Eqn  (21)  and  eliminating 
Fes  and  F^.^. ,  i.  e., 


L12B  Lji 


-4iL-2B-^a('P;2)-^U 


21 


-4iL;2A-'a('P-2)^‘L: 


21 


L^-Ll^A- 


-'21 


0 


-4,L;2A-'(l-/,a('P22)-') 


(29) 


Eqn  (29)  is  highly  nonlinear  because  the  matrices  Sq,A,  B,  'P22’  and  ^22 ,  are  all  nonlinear 
functions  of  the  crack  displacements  Ad;„j  (through  ao  and  ao  -  see  Eqn  26c  and  Eqn  28b),  which 
are  not  known  because  at  this  point  the  only  information  is  the  external  displacements  d~^,  and 
.  We  shall  defer  the  discussion  of  solution  procedure  until  section  3.3. 

3.2  A-FE  Formulation  with  Two  Interactive  Cracks 


The  internal  DoFs  associated  with  the  four  internal  nodes,  together  with  the  regular  nodal 
DoFs,  allows  for  stress  evaluations  at  each  subdomain,  making  it  possible  for  further  augmentation 
of  the  element  for  hosting  more  than  one  cohesive  crack.  A  secondary  crack  can  be  introduced  to 
a  subdomain  as  long  as  it  has  more  than  one  regular  (external)  node.  In  such  a  case,  there  are 
altogether  three  basic  cracking  configurations  as  shown  in  Fig  4a(l),  Fig  4a(2)  and  Fig4a(3),  where 
lei  is  the  first  crack  length  and  Ui  is  the  second  crack  length.  The  second  crack  intersects  and  cuts 
the  first  crack  into  two  segments  with  length  y/ei  (left)  and  (l-Y)/ei  (right).  The  three  subdomains 
are  noted  with  Q' and  Q.'  as  shown  in  Figure  4b. 

The  deformed  configurations  corresponding  to  the  three  cases  in  Figure  4(al-a3)  are  shown 
in  Figure  4(bl-b3).  To  describe  the  discontinuity  induced  by  the  second  crack,  four  more  internal 
nodes,  7,  8,  7’,  and  8’  are  introduced  and  a  local  coordinate  system  is  defined  by  the  angle  between 
its  tangential  direction  (S2)  and  global  x-direction,  02.  The  rotation  matrix  associated  with  this  crack 
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is  [R2]: 


cos  62  sin  62 
-sin^T  cos  6, 


With  the  8  internal  nodes,  it  is  possible  to  integrate  the  subdomain 


stiffness  matrix  for  each  subdomain.  In  the  following,  we  show  that  a  unified  elemental 
equilibrium  equation  can  be  derived  for  all  the  three  configurations. 


Figure  4  Illustration  of  a  2D  A-FE  with  two  intra-elemental  cracks  (3  subdomains) 


Despite  the  different  crack  configurations  in  Figure  4,  the  internal  DoF  and  force  arrays 
can  be  arranged  in  a  unified  way  as  below. 


■dsg.y,,  djjjj 


*786 


—  F  •  F''“  —  F  •  F'“  —  F 

*int  “^6'5'’  *int  “*58'7’’  *  int  “*786 


(30) 


This  enables  a  unified  expression  for  the  equilibrium  of  all  three  subdomains  as  follows 


(a) 

(b) 

(c) 

(d) 

(e) 

(f) 


(for  Q^) 
(forQ''  ) 
(forQ'  ) 


(31) 


Here,  the  sub-matrices  (L',.  ,Lf,  andLf  (/,  7  =  1,2))  are  the  stiffness  matrices  of  the 
subdomains  andQ^^  ,  respectively.  They  were  computed,  according  to  the  respective 
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crack  configurations  in  Fig  4(bl-b3),  with  the  sub-domain  integration  scheme  discussed  in  section 
3.1. 


To  integrate  the  cohesive  stress  along  the  three  cohesive  crack  segments  into  respective 
internal  nodal  forces,  the  first  order  Gaussian  Integration  scheme  with  one  integration  points  was 
used  on  each  of  the  three  cracked  segments  (i.e.,  assuming  uniform  cohesive  stresses  along  each 
crack  segment).  The  integration  points  are  at  the  center  of  each  crack  segment  as  indicated  in 
Figure  4(b)  by  the  solid  dots  labeled  by  I,  II,  and  III. 

Denote  the  cohesive  stresses  at  these  points  as  ( cr/'* ,  ),  ( ),  and  ( 

).  The  superscript  (i,  j,  k,  I,  m,  n)  are  the  free  indices  of  the  assumed  cohesive  segment  numbers. 
The  relation  between  these  cohesive  stresses  and  the  internal  nodal  displacements  is 

{r/'^ ,  }^  =  Ooj  +  ttoiRi  ((1  -  r  /  2)d6,  +  (r  /  2)d5. )  -  ao^R;  (d^  +  d^ )  /  2 

<  =  Oojj  +  ttojjRi  ((1  -  r  /  2)d6,  +  (1  +  r  /  2)d5, )  -  tto^R;  (d^  +  d^, ) /  2  (32a) 

=  <^0111  +  «omR2  (d^.  +  d3,)/2  -  ttojuR^  (d^  +  d3)/2 


where 


ttg,  =  Diag[«/‘V„*'^]; 

The  internal  forces,  {F^s. 

respective  crack  lengths  are 


'on 


ttoii  =  Diag[«/'V„*'’]- 


.},  and  {F^gg},  integrated 


^“ii  ~  J 

(32b) 

a,.„=Diag[«;'"V„'"^]; 

(32c) 

from  the  cohesive  stresses  along  the 


F  =  T^l 

^6'5'  ( 

‘■1  ’ 

(J) 

F  _  t'- 
^786  -■'R 

ir  (0 

^(7) 

(33a) 

F  -  T'^ 

^58'7'  -^R 

-{rr. 

^  ('w) 

where  ,  and  are  the  respective  integration  matrices  as  below: 


17 


r{2-r)h,K  . 

r  0  \(i-r)LK 

t^  =  t  rXK  42R2"  ;  Tr  =  -  (i-r)4.R/ 

0  0 


(33b) 


Substituting  Eqn  (32)  into  (33)  and  re-organizing  terms,  the  internal  forces  become: 


Sj  ttjj  aj2  «i3  dg,5 

~  “i  ^2  '  ®21  ®22  ®23  "  ^786 

So  tto,  a„  a„  d„, 


ttji  a32  tt33j  I  d58.7’ 


(34a) 


where 


S.  - 


,S2-T^r“ 

,*^011)  1*^0111 


,  S3^T- 


(34b) 


^  (2-r)aoiRi  ^ttoiRi 

“  2  (l-r)aonR,  (l  +  r)ao„R, 


U.  0  0  0 

®13  =  “Tjj 

2  -aoijRi  -ttoyRi  0 


1  d  ®qjRj  tXgjRj 

«?7  =  “  T„ 

2  _“®omR2  ~®omR2  d 

l_J(l-r)aonRi  (l  +  r)aonRi 

tto,  =  —  T„ 

^  2  L  0  0 

1  rpr-  “®0nRl  “®0IlRl  d 

«33=“T„ 

2  d  aQjjjR2  aQjjjR2 


2  "  0  0  0 


(2-^)aoiRi  xttoiRi 

a,i  =  —  T„ 

"‘2  0  0 


•  a  =— T 


®32  ~  2 


uJ«  » 


d  ®0I1iR2 

^0111^2  J 

0 

0 

.“®oinR2 

“®omR2 

(34c) 


Substituting  Eqn  (34a)  into  expressions  (b,  d,  and  f)  in  Eqns  (31)  and  solving  the  internal 
displacements  as  functions  of  the  regular  (external)  nodal  displacements. 


L22  (*11  (*12 

—(*21  L22  —  (*22 


S2  +  -L'l 


^22  ®33  [djg.y. 


The  internal  nodal  displacements  can  be  further  expressed  explicitly  as  functions  of 
external  nodal  displacements  as  below: 
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(36) 


% 

Ys 

y; 

dvse 

>  = 

Y4 

Y5 

Y, 

^dss'V’, 

Y, 

Y2 

Y3_ 

1 

X,; 

X„ 

X.2 

< 

dZ 

J 

00 

1 _ 

d:. 

where  the  matrices  (p  =  1,  2, ..  15)  and  =  1,  2, ..  9)  can  be  found  in  [21]. 

Finally,  by  substituting  Eqn  (36)  into  Eqn  (31a,  31c,  and  31e)  to  eliminate  the  internal 
nodal  displacements,  the  fully  condensed  elemental  equilibrium  equation  is  obtained  as 


^11  ^12X10 

i^r2x„ 

r  3 

dL 

F^“ 

ext 

L^Ys 

l^y; 

L^Xv 

L^+L-X, 

L-X, 

< 

dZ 

.  =  - 

ext 

>  — 

Lr2Y, 

Lr2Y3 

< 

S2 

L^Xn 

L;2Xi4 

l;i+l;2x,5 

d:. 

F^ 

ext 

L^Yv 

L^Ys 

L^Y, 

S3, 

(37) 


Note  that  Eqn  (36)  and  Eqn  (37)  share  the  same  matrices  (p  =  1,  2,  . .  .15),  (  ^  =  1, 

2, ...9),  ),  and  Sr  (r  =  1,  2,  3),  which  are  all  functions  of  the  characteristic  cohesive  stresses 
(Ooj,Oq;j,  and  (Eqn  32b)  and  the  cohesive  stiffnesses  (aopa^jp  and 

(Eqn  32c).  They  are  not  known  at  this  point  because  the  crack  displacements  are  not  known  yet. 
The  solving  procedure  for  these  equations  will  be  discussed  in  next  section  3.3. 


3.3  Consistency-check  Based  Elemental  Condensation  Algorithm 


In  section  3.1  and  3.2,  the  elemental  equilibrium  equations  for  single  crack  (Eqn  29)  and 
for  two  interactive  cracks  (Eqn  37)  have  been  derived.  However,  these  nonlinear  equations  have 
yet  to  be  solved.  In  literature,  Eqn  (29)  or  Eqn  (37)  is  typically  re-formulated  into  an  incremental 
form  and  then  solved  with  either  the  direct  matrix  inversion  technique  [50,  51,  53]  or  a  Newton- 
Raphson  type  of  iterative  solver  [49,  75,  76].  This  is  known  as  the  static  condensation  process. 

In  this  paper,  the  explicit  nature  of  the  solutions  enables  a  very  efficient  elemental 
condensation  process  that  is  completely  different  from  what  have  been  used  in  literature.  Instead 
of  solving  the  incremental  form  of  the  elemental  equilibrium  of  Eqn  (29)  or  Eqn  (37),  in  this  paper 
we  first  solve  the  internal  nodal  displacements  (Eqn  (28)  or  Eqn  (36))  analytically  through  a  simple 
consistency  checking  procedure.  It  will  be  shown  that  once  these  equations  are  solved,  the 
equilibrium  equations  (Eqn  (29)  or  (37))  are  satisfied  with  mathematical  exactness  by  substituting 
the  related  matrices  into  them. 


Single  Intra-Elemental  Crack 

In  this  case,  the  internal  crack  displacements  have  been  given  in  Eqn  (28).  It  is  noted  that, 
once  Eqn  (28a)  is  solved,  i.e.,  So,A,  B,  'P22’  established,  simply  substituting  these 

matrices  into  Eqn  (29)  will  guarantee  the  elemental  equilibrium.  We  first  convert  the  global  crack 
displacements  in  Eqn  (28a)  into  local  crack  displacements  using  Eqn  (26b),  which  gives 
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^+(b-'  -4A^V^22)-')L2,d:,, +(4,B-y'p;,)-'  -  A-')L;,d:,,^ 

(38) 

In  Eqn  (38)  the  superscripts  {i,i,  k,  1)  in  the  right-hand- side  are  included  to  emphasize  that 
the  matrices  Sq,A,  B,  ^22,  and'P22  ^  functions  of  the  cohesive  slope  matrix  and 

characteristic  stress  matrix  (Eqn  26c  &  28b).  Thus,  for  any  trial  cohesive  segment 

combination  (/,  j,  k,  1),  a  set  of  local  crack  displacements  is  immediately  solved  from  Eqn  (38). 
The  solved  local  crack  displacements  are  then  subjected  to  a  consistency  check,  and  the  true 
solution  is  found  if  the  following  statement  is  true 

“  (0  —  (7)  —  (k)  —  (/) 

e  5s  is  true)  &&  ((5’„66.  e  5n  is  true)  &&  {5^^^  e  5s  is  true)  &&  {5^^^  e  5„  is  true) 

~  (0  —  (7)  —  (*:)  —  (0 

where  5s  ,  5n  ,  5s  and  5s  are  the  respective  crack  displacement  ranges  of  the  assumed 
segments  (ij,  k,  1)  (Section  2.2). 

This  is  the  core  idea  behind  the  solving  algorithm  used  in  this  paper.  The  mathematical 
exactness  (in  piece-wise  linear  sense)  of  the  solution  is  guaranteed  by  the  solution  uniqueness  of 
the  local  problem  (the  solution  uniqueness  for  linear  material  with  nonlinear  discontinuities  has 
been  discussed  in  [77]).  The  detailed  algorithm  formulation  and  extended  discussions  on  its 
numerical  performance  in  solving  fracture  problems  can  be  found  in  [78]. 

Two  Intra-Elemental  Cracks 


'^s66' 

^„66' 


^555' 


=  Diag[R/;R/] 


In  this  case,  the  crack  displacements  as  explicit  functions  of  external  displacements  (such 
as  the  Eqn  (38)  for  the  single  crack  case)  are  not  readily  available.  However,  the  solving  procedure 
presented  above  remains  applicable  in  solving  the  internal  displacement  equation  of  Eqn  (36a). 
The  equation  is  repeated  below  with  superscripts  added  to  emphasize  the  dependence  of  (p  = 

1,2,..  15),  Y^{q=l,2, ..  9),  and5r(r=  1,  ..3)  on  the  cohesive  segment  combination  {i,  j,k,l,m,n) 


de'S' 

f 

1 

Ys 

1 

'Si' 

A 

Xi3 

Xi4 

Xl5" 

d786 

>  = 

Y4 

Y5 

Y<, 

« 

S2 

> 

+ 

X,o 

X„ 

X,2 

< 

d:;, 

dsS'V’, 

V 

Y^ 

Y3_ 

S3. 

J 

1 

-J 

Xs 

x,_ 

d;, 

(39) 


Eor  a  trial  cohesive  segment  combination  {i,  j,k,l,m,n)  and  given  external  nodal  displacements, 
a  set  of  solution  to  the  internal  nodal  displacements,  ‘le'S’  “  ’^786  -  {d7’dg,dg}  ,  and 

djg.y,  =  {d5,dg,,d2,}  ,  is  immediately  solved  from  Eqn  (39).  Then,  the  crack  displacements  at  the 
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integration  point  of  each  segment  (points  I,  II,  and  III  in  Figure  4)  are  computed  and  converted  to 
local  coordinates,  i.  e.. 


'  p,, ,  J"  =  r;  (a  -  /  /  2)d,  +  (r  /  2)d,,  -  (d, + d, )/ 

.  f  =  R,"  (d  -  y / 2)d,,  +  (1  +  y) / 2d,,  -  (d,  +  d,, ) / 2)"  ''‘  '''""’  (40) 

{^,m.  <5.m  r  =  R/  ((d,,  +  d,. )  /  2  -  (d,  +  d.)  /  2f 

The  local  crack  displacements  are  then  subjected  to  a  consistency  check  similar  to  the  single  crack 
case  and  the  exact  solution  is  found  when  the  following  statement  is  true. 


e  S.!'  ^  true)  &&  (S,^j  e  true)  &&  (S^„  e  Ss^  ^  true)  &&  e  S,!  *  true) 

&&  e  s!'  *  true)  &&  e  S,!  ^  true) 

Once  the  solution  is  found,  (p  =  1,  2,  ...15),  =  1,  2,  ...9),  and5r(r  =  1,  2,  3)  are 

all  established  with  the  determined  (i,  j,  k,  I,  m,  n).  Substituting  them  into  Eqn  (29)  the  fully 
condensed  elemental  equilibrium  is  then  satisfied  with  mathematical  exactness  (in  piece-wise 
linear  sense  and  within  the  limit  of  FE  discretization  accuracy). 

3,4  Numerical  Implementation 

The  new  A-EEM  and  the  solving  procedure  above  can  be  integrated  into  any  standard  FE 
programs  as  a  general  purposed  element  because  the  elemental  locality  is  completely  retained. 
Here  we  briefly  describe  how  to  implement  it  into  an  existing  EE  program  as  an  element. 

Eor  any  displacement  based  EEM  programs,  the  external  nodal  displacements  are  passed 
into  the  element  to  establish:  1)  instantaneous  element  stiffness  matrix  or  Jacobian  matrix,  and  2) 
the  external  force  array  (also  called  right-hand-side  or  RHS).  Both  matrices  have  been  derived  and 
solved  explicitly  in  section  3.1  for  single  crack  case  (Eqn  29)  and  in  section  3.2  for  two  intra- 
elemental  cracks  (Eqn  37). 

The  local  crack  displacements  were  saved  as  solution-dependent  variables  (or  state 
variables)  and  subjected  to  constant  updates  so  that  their  maximum  values  ever  experienced  were 
always  available  for  assessing  loading  or  unloading  condition  (for  consistency  check).  Another 
minor  issue  is  to  guarantee  the  deformation  continuity  across  the  element  edge  that  is  shared  by  a 
crack-tip  element  and  the  element  immediately  ahead  of  it,  which  is  a  standard  continuous  element. 
This  was  conveniently  done  in  this  paper  by  enforcing  the  internal  nodal  displacements  at  the  edge 
to  be  dj  =  dj,  =  N ,  where  5  and  5’  are  the  internal  nodes  at  the  crack-tip;  A  and  B  are 

the  two  external  nodes  of  the  element  edge  hosting  the  crack-tip;  and  Na  and  Nb  are  the  standard 
line  shape  functions  associated  with  node  A  and  B. 

One  of  the  advantages  of  this  method  is  that  it  does  not  need  a  sophisticated  global  tracking 
algorithm  such  as  the  level-set  functions  in  the  X-EEM  for  a  propagating  crack,  which  requires 
heavy  calculations  when  implemented  in  a  narrow  band  manner.  However,  it  is  beneficial  to  record 
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the  crack  path  so  that  the  inter-element  continuity  of  a  crack  is  reinforced.  This  was  done  with  a 
common  block  accessible  to  all  user  elements,  within  which  the  geometric  information  of  a  crack 
is  constantly  updated  as  it  propagates.  A  detailed  data  structure  and  propagation  rules  regarding 
multiple  crack  interaction  has  been  given  in  detail  in  [21]. 

We  note  that  the  new  A-FEM  only  needs  to  record  the  cracks  developed  in  response  to  the 
local  stress  conditions.  There  is  no  need  to  actively  assign  an  incremental  crack  length  for  crack 
extension.  The  crack  initiation  in  an  element  is  purely  determined  by  the  local  stresses  in  the  crack 
front  element  accordingly  to  a  specified  initiation  criterion  (such  as  the  maximum  principal  stress 
criterion  used  in  the  following  numerical  examples  in  section  4).  The  propagation  is  purely 
determined  by  the  cohesive  crack  propagation  rule  of  Eqn  (15). 

4.  3-D  A-FEM  FORMULATION 

Extension  of  the  2D  A-EEM  formulation  to  3D  is  relatively  straightforward.  Here  we 
briefly  present  the  3D  A-EEM  for  the  4-node  tetrahedron  element  as  shown  in  Eigure  5  (a).  [79]. 
There  are  two  possible  configurations  if  the  element  is  cut  by  a  crack  surface;  i)  the  element  is  cut 
into  one  tetrahedron  subdomain  and  one  wedge  subdomain  as  shown  in  Eig  5(b),  and  ii)  the 
element  is  cut  into  two  wedge  subdomains  as  shown  in  Eig  (5c). 

The  tetrahedron- wedge  configuration  with  the  triangular  crack  plane  (Eigure  5(b))  has  6 
internal  nods  as  5,  6,  7,  5’,  6’  and  7’  while  the  wedge-wedge  configuration  with  quadrilateral  crack 
plane  has  8  internal  nodes  5,  6,  7,  8,  5’,  6’, 7’  and  8’.  Displacement  jumps  across  the  cohesive  crack 
is  simply  the  difference  in  displacements  between  node-pair  5-5’,  6-6’,  7-7’,  and  8-8’.  Defining  a 
local  coordinate  system  as  shown  in  Eig  5(d)  with  the  surface  out-normal  (N  )  as  the  local  z-axis, 
the  rotation  matrix  can  be  obtained  as 


zrZj 

n, 

R  = 

h 

m2 

(41a) 

h 

m3 

where  /,,m  .  and  (i,y,  k  =  1,2,3)  are  the  direction  cosines  of  the  crack  plane  with 

N  =  /j/  +mj  -\-n^k  being  the  out-normal  direction,  T  =  1^1  +^1^]  +nj<.  and  S  =  njZ  -\-n2j  +  n^k 
being  the  two  orthogonal  in-plane  shear  directions,  and 

/,  =  Cosa.,  m.  =  Cosp.,  n.  =  Cosy.  (/  =  1,2,3) .  (41b) 

a.,  p.,  and  are  the  angles  between  the  local  axes  and  global  axes  as  illustrated  in  Eigure  5(d). 

Eollowing  the  similar  procedure  of  2D-AEEM  in  section  3.1,  the  local  equilibrium  equation 
can  be  written  as  the  identical  form  of  Eqn  (19): 
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The  only  differences  between  Eqn  (19)  and  (42)  are  the  dimension  of  the  submatrices  and  the 
definition  of  external  and  internal  DoFs  and  forces.  For  the  tetrahedron-wedge  cut  configuration 
(Fig  5  a), 

-F-  F  -F  •  -H  •  H  -H 

*ext  *4’  *ext  *123’  **ext  **4  ’  **ext  **123 

F^  -F  •  F  -F  •  -H  •  d  -H 

*int  *5'6'7'’  *int  *567’  **int  **5'6'7' ’  **int  **567 

and  for  the  wedge- wedge  cut  configuration  (Fig  5b), 

F^  -F  •  F  -F  •  -d  •  d  -d 

*ext  *34’  *ext  *12’  **ext  **34’  **ext  **12 

F^  -F  •  F  -F  •  d^  -d  d  -d 

*int  *5'6'7'8'’  *iiit  *5678’  **iiit  **5'6'7'8' ’  **int  **5678 

The  internal  force  array  integrated  from  cohesive  stresses  between  the  cohesive  crack 
surfaces  can  be  written  as 

F.;  =  -K  =  A  (So  +« •  )  ,  (43a) 


where  Ae  is  the  area  of  the  crack  surface.  S„  and  a  are  the  characteristic  stress  and  cohesive 
stiffness  matrix  similar  to  those  in  Eqn  (26).  For  the  tetrahedron- wedge  cut  configuration. 


S„  =Diag[R^;R^;R^](T,„J^a„ 
a  =  Diag[R^  R^;  R^]  (T^„J^  Diag[R;  R;R] 


tto  =  Diag[a® ;  a® ;  a® ;  af * ;  ;  Af  af ;  a® ;  a®  ] 


(43b) 


where  and  are  the  respective  integration  and  interpolation  matrix  for  the  3-node  cohesive 

crack  surface.  Each  of  the  free  indices  (iJ,k,l,m,n,o,p,q)  can  range  from  1-4  for  the  cohesive  law 
shown  in  Figure  2. 

For  the  wedge-wedge  cut  configuration. 


So  =Diag[R®R®R";R"](T,„j"tT„ 
a  =  Diag[R^  R^;  R^;R^]  (T^„J^  ao  Diag[R;  R;R;R] 


- 


On  =  1  r;  -  a 


(0^(0.  .f(;)  _/yO')^(y).^(^) 

’>^t  ^  ’ 


tto  =  Diag[af ;  a® ;  a® ;  af ;  ^ ] 


(43c) 


In  this  case  the  respective  integration  and  interpolation  matrix  for  the  4-node 

cohesive  surface.  Each  of  the  free  indices  {ij,k,l,m,n,o,p,q,r,s,t)  can  range  from  1-4  for  the  cohesive 
law  shown  in  Figure  2. 
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With  the  characteristic  stress,  cohesive  stiffness,  interpolation  and  integration  matrices 
fully  defined,  the  equation  that  relates  internal  displacements  to  the  external  displacements,  and 
the  equation  of  elemental  equilibrium  are  found  to  be  in  exactly  the  same  matrix  form  of  Eqn  (28) 
and  (29),  respectively.  Thus  the  consistency-check  based  procedure  for  elemental  condensation  as 
described  in  section  3.3  can  be  readily  applied  to  solve  the  element  equibrium. 

One  of  the  biggest  challenges  in  3D  fracture  modeling  is  to  maintain  the  conformity  of  the 
evolving  crack  surface.  A  crack  surface  tracking  scheme  is  absolutely  necessary  and  proper 
numerical  treatments  have  to  be  done  to  guarantee  that  the  crack  surface,  which  may  cut  multiple 
elements  before  they  merge  at  some  point,  maintains  a  relatively  smooth  planar  configuration 
without  too  much  out-of-plane  oscillation.  Several  crack  tracking  schemes  such  as  local  and  global 
tracking  algorithms  have  been  developed  in  the  literature  [80-83].  Here  we  adopt  the  local  tracking 
algorithm  first  proposed  by  Aries  and  Belytschko  [80].  The  detailed  algorithm  will  not  be  given 
here  but  interested  readers  are  referred  to  [80] . 


Figure  5:  Element  illustration  for  (a)  a  regular  4-node  tetrahedron  element  with  two  possible  cut  configurations 
by  a  crack  plane,  (b)  the  A-FE  with  tetrahedron-wedge  configuration,  (c)  the  A-FE  with  wedge-wedge 
configuration,  (d)  definition  of  local  coordinates  for  the  crack  plane. 


5.  DEALING  WITH  MATERIAL  ORTHOTROPY  AND  NONLINEARITY 

A  majority  of  modern  composites  are  highly  orthotropic  due  to  the  directional 
reinforcement  structures.  The  material  orthotropy  is  not  only  reflected  in  the  orthotropic  elastic 
constants,  but  also  in  their  progressive  damage  modes,  i.e.,  matrix  tension  and  compression  failure, 
fiber  direction  tensile  rupture  and  compressive  kinking,  shear  failure,  and  interlaminar 
delamination.  The  various  failure  modes  can  appear  simultaneously  but  their  interactions  have  not 
been  fully  studied  due  to  lack  of  analytical  tools.  Furthermore,  it  has  been  shown  that  typical 
polymer  matrix  composites  exhibit  various  nonlinearity  at  continuum  level  (i.e.,  before  initiation 
of  discrete  cracks).  These  include  strong  shear  nonlinearity,  asymmetric  tensile  and  compressive 
modulus,  and  often  times  with  strong  nonlinear  stress-strain  relation  under  longitudinal  and 
transverse  compression  [29,  84,  85].  While  the  strong  influence  of  the  shear  nonlinearity  on 
interlaminar  delamination  has  been  well  appreciated  [9,  13,  20,  86,  87],  overall  the  complex 
interaction  between  such  material  nonlinearity  and  the  intra-ply  damage  initiation  and  propagation 
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has  not  been  well  studied,  largely  due  to  the  lack  of  high-fidelity  analytical  means  to  couple  them 
into  a  single  analysis.  In  this  study,  the  following  nonlinearities  were  explicit  considered: 

5.1  Asymmetric  Tension/Compression  in  Longitudinal  and  Transverse  Direction 

Based  on  experimental  evidence,  the  stress-strain  relation  under  longitudinal  or  transverse 
tension  is  typically  linear  elastic  until  the  respective  strength  value  (At  or  Ft)  is  reached.  This  is 
represented  with  a  tensile  modulus,  i.e.,  for  longitudinal  tension  or  for  transverse 

tension.  Under  longitudinal  or  transverse  compression,  however,  typical  PMCs  exhibit  relatively 
strong  nonlinearity  in  the  stress-strain  relation  up  to  the  respective  strength  value  (Ac  or  Fc),  as 
shown  in  Figure  6(a).  In  this  paper,  the  full  stress-strain  relation  is  written  as  follows 

\^kT  0 

si  <s^<Q  {k  =  l  or  2)  ^-44^ 

^k  <  4 


^k  {^k )  = 


^kC  ^k 


^kC^k 


^  WkC^k  !  ^k 


0 


where  the  subscript  k  =  1  stands  for  longitudinal  response  and  k  =  2  refers  to  transverse  response. 
Ekc  is  the  compressive  modulus  and  si  is  the  elastic  limiting  strain  in  compression.  (7°  is  an 

asymptotic  stress  and  \  is  a  power  index.  They  control  the  shape  of  the  nonlinear  constitutive 
relation  and  have  to  be  characterized  from  experiments.  The  nonlinear  constitutive  relation  was 
implemented  as  a  damage  model  (Chapery  model),  i.e.,  the  instantaneous  secant  modulus  derived 
from  Eqn  (44)  was  used  to  form  the  material  stiffness  matric  (D  matrix  in  Eqn  18)  and,  upon 
unloading  at  a  strain  si  >  si ,  the  unloading  path  points  directly  towards  the  origin  (Eig  6a). 


5.2  Shear  Nonlinearity 

Similar  nonlinear  model  was  used  to  account  for  the  shear  nonlinearity,  as  shown  in  Eigure 
6(b).  The  shear  stress-strain  curve  can  be  written  as: 


'-12 


(^2)  = 


Gi2  F2 


^2 


1  +  RF2/^I2 


0  "s 


|F2|^n2 

1^2 1  ^  Tn 


(45) 


Where  G°2  is  the  elastic  shear  modulus  and  /l^  the  elastic  limiting  shear  strain.  is  an 
asymptotic  shear  stress  and  is  a  power  index.  They  control  the  shape  of  the  nonlinear 

constitutive  relation  and  have  to  be  fitted  from  experimental  data.  The  nonlinear  shear  stress-strain 
relation  was  also  implemented  as  a  damage  model  (Chapery  model),  i.e.,  the  instantaneous  secant 
modulus  derived  from  Eqn  (45)  was  used  to  form  the  material  stiffness  matrix  (i.e.,  the  D  matrix 
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in  Eqn  18),  and  upon  unloading  at  a  strain  >  the  unloading  path  points  directly  towards 

the  origin  (Fig  6b). 


<7* 

X,(orY,) 

e; 

^kT 

1  ✓  1 

E  ^ 

^kC 

-X^{or  -FJ 

(a) 

Figure  6  (a)  asymmetric  tension/compression  stress-strain  relation,  and  (2)  shear  nonlinearity. 


5.3  Cohesive  Laws  for  Fiber  Rupture,  Fiber  Kinking,  and  Matrix  Cracks 

Similar  to  all  other  numerical  techniques  such  as  X-FEM  and  PNM-FEM  that  uses  cohesive 
zone  models  to  describe  the  progressive  failure  processes,  in  this  study  it  is  assumed  that  once  an 
initiation  criterion  of  a  certain  damage  mode  is  satisfied,  a  cohesive  crack  with  corresponding  type 
is  inserted.  Further  cohesive  damage  evolution  until  the  crack  becomes  traction  free  is  based  on 
the  normal  and  shear  traction-separation  laws  and  a  mixed  mode  propagation  criterion.  As 
discussed  earlier,  it  is  very  critical  to  consider  tensile  and  compressive  cohesive  responses 
explicitly  and  simultaneously,  and  it  is  also  very  important  to  consider  the  possible  coupling 
between  the  two  modes.  In  this  study,  we  propose  the  unified  tension-compression  cohesive  laws 
for  normal  cohesive  response  and  apply  a  commonly  used  shear  cohesive  law  as  shown  in  Figure 
7.  Triangular  shaped  cohesive  laws  were  chosen  for  all  modes  for  their  simplicity. 

The  shear  (mode  II)  response  shown  in  Figure  7(a)  is  a  standard  one  used  widely  in 
literature  (e.g.  [67,  69]).  The  response  is  anti-symmetrical  with  respect  to  positive  and  negative 
sense.  The  two  governing  cohesive  parameters  are  the  shear  strength  5'j2  and  a  mode  II  toughness 

Tj ,  from  which  the  critical  shear  displacement  can  be  obtained  as  /  3^2 . 

The  fiber  direction  (or  longitudinal)  tensile  and  compressive  cohesive  responses  are  unified 
into  a  single  cohesive  law  as  shown  in  Figure  7(b)  for  tensile  rupture  and  compressive  kinking. 
The  tensile  response  is  governed  by  the  longitudinal  tensile  strength  and  a  toughness  T^j, , 

from  which  the  critical  displacement  can  be  obtained  as  -  lYp^  /  ;  while  the  compressive 

response  is  governed  by  the  longitudinal  compressive  strength  -X^  and  a  toughness  T^^  from 
which  the  critical  compressive  displacement  can  be  obtained  as  -  -2Ypp  /  X^ . 
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Similarly,  the  cohesive  model  for  matrix  tension  and  compressive  responses  are  unified 
into  a  single  cohesive  law  as  shown  in  Figure  7(c).  The  tensile  response  is  governed  by  the 
transverse  tensile  strength  and  a  toughness  F^j, ,  from  which  the  critical  displacement  can  be 

obtained  as  =  2r„j,  /  .  The  compressive  response  is  governed  by  the  transverse  compressive 

strength  -Y^  and  a  toughness  F^^  from  which  the  critical  compressive  displacement  can  be 
obtained  as 


Figure  7  cohesive  laws  for  (a)  (matrix)  shear  damage  mode;  (b)  fiber  tensile  rupture/compressive  kinking  damage 
mode;  c)  matrix  tension/compression  damage  mode.  Here  the  first  subseript  “F”  denotes  variables  assoeiated 
with  fiber  rupture  or  kinking  damage,  “M”  refers  to  variables  assoeiated  with  matrix  tension  or  eompression 
damage.  The  seeond  subseript  “T”  or  “C”  differentiate  the  tensile  (“T”)  and  eompressive  (“C”)  part  of  the 
cohesive  laws. 

We  note  that  the  negative  normal  displacements  (^mc’^fc)  the  compressive  cohesive 
laws  are  not  unphysical  for  typical  composite  failure  under  compression  dominant  conditions:  the 
compressive  kinking  and  compressive  matrix  crushing  do  induce  finite  amount  of  local  negative 
normal  displacements,  if  the  associated  damage  bands  are  modelled  with  cohesive  cracks  with 
zero-thickness  cohesive  zones,  as  implemented  in  this  and  other  studies.  Numerically,  such 
negative  normal  displacements  lead  to  mesh  interpenetration  if  they  are  modelled  as  explicit 
cohesive  elements,  causing  numerical  divergence  (negative  Jacobian).  However,  this  is  not  an 
issue  with  our  embedded  cohesive  zone  model  in  the  A-FEM. 

It  is  noted  that  in  Figure  7,  all  the  cohesive  laws  are  piece-wise  linearized.  Each  of  the 
segment  with  a  constant  slope  is  indexed  with  a  number.  This  is  done  in  accordance  with  the 
consistency-check  based  condensation  algorithm  discussed  in  section  3.3.  Each  segment  is 
characterized  by  a  critical  stress  ( f  or  ),  a  corresponding  local  crack  displacement  ( 

or  ),  and  a  constant  stiffness  (a/'^or  cCpiJ'^'*)-  Here  the  superscripts  i  (=1,2,3)  andj  (=  1,  2, 

3;  r,  2’,  3’)  are  free  indices  for  the  cohesive  segments  with  constant  stiffness. 

For  each  cohesive  law  in  Figure  7,  segment  3  and  6  are  defined  for  irreversibility.  The 
characteristic  stresses  and  displacements,  which  are  marked  with  superscript  are  the  maximum 
historical  values  ever  experienced  during  the  entire  deformation  history  up  to  this  point.  Should 

unloading  occur,  it  follows  an  instantaneous  path  with  a  constant  stiffness  of  /  ^rim*  ■  this 
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study,  these  historical  quantities  are  recorded  for  each  damage  modes  and  they  play  key  roles  in 
enabling  partial  damage  coupling  to  be  discussed  shortly. 

5.4  3D  Geometric  Models  for  Binary  Model  Analysis  of  Textile  Composites 

To  handle  large  textile  CMC  structures,  we  have,  follow  the  same  augmentation  process  as 
in  section  2,  developed  a  breakable  ID  tow  element  (Figure  8),  which  can  be  coupled  with  the  3D- 
A-FEM  to  simulate  progressive  damage  and  the  failure  of  reinforcement  tows.  The  ID  tow 
elements,  which  follow  the  spatial  loci  of  each  of  the  tow  center  lines,  are  embedded  in  a  3D 
effective  medium,  which  is  modelled  by  3D  A-FEM  elements,  a  strategy  that  is  based  on  the  Binary 
Model  (BM)  formulation  established  in  prior  work  [88-90].  Possible  tensile  tow  rupture  and 
compressive  tow  kinking  are  incorporated  via  a  ID  cohesive  law. 


(a)  (b) 

Figure  8.  (a)  a  simple  BM  model  containing  a  prism  EM  and  a  straight  tow  (red  dashed  line),  (b)  The  simulated 
elemental  stress-strain  response  of  the  BM  (red-solid  line),  which  is  a  summation  of  the  solid  element  responses 
(blue  dashed  line)  and  the  tow  response  (green  dashed  line). 

We  have  completed  and  improved  the  3D  A-FEM  with  shear  nonlinearity  to  simulate 
realistic  3D  textile  composites  under  the  BM  scheme.  The  subcontractor  of  this  project.  Dr.  Brian 
Cox  has  developed  a  very  efficient  algorithm  for  generating  tow  elements  with  realistic  spatial 
distributions  in  textile  composites  [91].  The  generator  converts  machine  instructions  for  a  3D 
weaving  loom  into  topological  ordering  rules,  which  are  then  converted  into  3D  spatial  models  of 
the  tows  in  the  textile.  The  input  deck  of  the  generator  consists  entirely  of  arrays  of  bounded 
integers,  which  are  defined  to  correspond  to  weaving  machine  actions,  thereby  guaranteeing  that 
any  choice  of  the  integers  is  in  principle  an  achievable  architecture  on  the  machine.  The  fact  that 
only  bounded  integers  are  used  as  input  allows  the  generator  to  be  linked  easily  with  an 
optimization  search  engine,  e.g.,  based  on  a  genetic  algorithm  or  an  ant  colony  search  algorithm 
Such  an  example  is  given  in  Figure  9. 

By  embedding  such  tow  elements  into  the  3D  effective  medium  modelled  with  our  3D  A- 
FEs,  very  encouraging  results  can  be  obtained.  One  of  such  example  will  be  detailed  in  next 
section. 
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Figure  9.  Illustration  of  the  tow-mesh  generator  the  geometric  model  generator,  (a)  Machine  actions  on  a  3D  weaving 
loom,  which  are  represented  as  arrays  of  bounded  integers  (b),  from  which  a  complex  3D  component  is 
automatically  generated  (c).  Nothing  in  the  input  deck  tells  the  generator  explicitly  that  it  is  making  a 
sandwich  structure. 


6.  NUMERICAL  EXAMPLES  AND  DISCUSSIONS 
6.1  Four-Point  Shear  Beam  Test 

In  this  section,  we  test  the  new  A-FEM’s  capability  in  simulating  the  crack  propagation  in 
a  4-point  shear  beam  test  reported  by  [92] .  The  problem  has  been  simulated  successfully  using 
extended  finite  element  method  (X-FEM)  by  Moes  and  Belytschko  [66].  Here  we  adopt  the 
numerical  model  setup  used  in  [66],  which  is  reproduced  in  Figure  10.  The  geometry  dimensions 
are:  b  =  200  mm;  I  =  4/?;  a  =  02b\  cib  =  0.4;  h  =  h  =  20  mm;  t  =  100  mm. 

Maximum  principal  stress  criterion  was  used  for  crack  initiation.  That  is,  crack  initiation 
from  a  pristine  elastic  element  occurs  when  the  maximum  principal  stress  averaged  within  the 
element  ) ,  reaches  the  mode  I  cohesive  strength  (cr) .  The  crack  path  is  perpendicular  to  the 
maximum  principal  stress  direction.  The  concrete  material  properties  are:  E  =  28000  N/mm^;. 
V  =  0.1;  Tj  =  0.145  N/mm;  <7  =  2.4  N/mm^ .  The  cohesive  law  used  is  triangular  type  with  an 
initial  slope  of  =  2.0  x  10"^  N/mm^  and  a  softening  slope  of  =  ~20  N/mm^ . 
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As  in  [66],  the  maximum  principal  stress  criterion  was  used  for  crack  initiation.  That  is, 
crack  initiation  from  a  pristine  elastic  element  occurs  when  the  maximum  principal  stress  averaged 
within  the  element  ,  reaches  the  mode  I  cohesive  strength  (cf) .  The  crack  direction  is 

perpendicular  to  the  maximum  principal  stress  direction. 

To  investigate  the  mesh  sensitivity  of  the  new  A-FEM,  the  problem  was  analyzed  by  five 
different  structured  meshes  with  characteristic  mesh  size  of  =  2,  4,  8,  13,  and  20  mm  (the  total 
elements  are  10760,  2849,  1094,  640,  and  410,  in  the  same  order),  and  three  unstructured  meshes 
with  /z  ~  4,  8,  and  15  mm  (total  elements  4062, 1333,  and  463,  respectively).  The  two  finer  meshes, 
h  =  2  and  4  mm,  are  typical  of  those  used  in  literature  [51,  52,  66,  75],  while  the  larger  meshes  (8, 
13.3,  and  20  mm)  are  used  in  this  study  to  explore  the  mesh  limit  of  the  new  A-FEM.  The 
unstructured  meshes  of  ~  4,  8,  and  13  mm  were  intended  to  check  the  mesh  sensitivity  and  the 
robustness  of  the  new  A-FEM.  Five  of  the  discretized  models  (structured  h  =  2,  S,  and  20  mm 
meshes  and  unstructured  h  ~  S  and  13  mm  meshes)  under  deformed  states  are  shown  in  Figure  10 
(b,  c,  d,  e,  f)  with  the  crack  trajectories  roughly  following  the  center  lines  of  the  white  bands. 

The  loading  parameters  were  set  to  be  identical  for  all  the  cases,  with  a  maximum 
prescribed  displacement  of  1.0  mm,  a  suggested  initial  incremental  size  of  0.001  mm,  and  a 
maximum  incremental  size  of  0.01  mm.  The  actual  size  of  each  increment  was  determined  by  the 
automatic  incrementation  scheme  in  ABAQUS  which  depends  on  the  convergence  rate. 


All  simulations  were  run  on  a  Dell  precision  M4600  (x64  bit)  mobile  workstation  with  Intel 
Core  17-2860  QM  CPU  @  2.5  GHz  and  with  8  GB  of  RAM.  Simulations  were  terminated  when 
the  load-point  displacement  reached  0.1  mm  (after  snap-back).  RIKS  option  (arc-length  method) 
was  invoked  in  attempt  to  capture  the  snap-back  behavior. 


(b)  2  mm  -  structured  (10760  elements) 


(c)  8  mm  -  structured  (1094  elements) 


(d)  20  mm  -  structured  (410  elements) 


(e)  8  mm  -  unstructured  (1233  elements) 


(f)  13  mm  -  unstructured  (463  elements) 


Figure  10  (a)  Numerical  models  and  simulated  crack  paths  for  structured  2mm  (b),  8  mm  mesh  (c),  and  20  mm  mesh 
(d),  and  unstructured  8  mm  (e)  and  13  mm  (f)  meshes. 


Numerical  Accuracy  and  Mesh  Sensitivity 
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The  simulated  load-displacement  curves  are  plotted  in  Figure  1 1(a).  The  X-FEM  prediction 
obtained  by  Moes  and  Belytschko  [66],  which  used  a  triangular  mesh  of  h  ~  3  mm,  was 
superimposed  on  this  plot  for  comparison.  For  each  case,  the  load  linearly  increases  with  the 
applied  displacement  until  crack  initiation  from  the  pre-crack  tip,  which  occurs  at  an  applied 
displacement  of  ~  0.04  mm.  After  the  initiation  the  load-displacement  curve  becomes  increasingly 
nonlinear  due  to  continuous  crack  growth.  The  peak  load  is  reached  at  the  displacement  of  about 
-0.09  mm.  After  the  peak  load,  all  curves  exhibit  a  strong  snap-back  behavior.  In  this  study,  the 
strong  snap-back  behavior  was  captured  using  the  arc-length  method  available  in  ABAQUS 
(Option  RIKS). 

The  A-FEM  computed  load-displacement  results  are  very  consistent  with  the  X-EEM  result 
obtained  in  [66].  In  particular,  the  load-displacement  curves  obtained  with  smaller  mesh  sizes  (i.e., 
h  =  2  mm,  4  mm,  and  8  mm)  are  almost  identical  with  the  X-FEM  results.  Eor  the  larger  meshes, 
i.e.,  h=  13  mm  and  20  mm,  the  curves  deviate  mildly  from  the  finer  mesh  curves  due  to  increases 
in  initial  stiffness.  The  systematic  increase  of  initial  stiffness  with  the  increase  of  mesh  sizes  is 
largely  due  to  the  inherent  numerical  inaccuracy  associated  with  the  4-node  plane  elements  (it  is 
well-known  that  larger  EE  elements  lead  to  overestimate  of  elemental  stiffness).  If  the  deviation 
of  initial  slope  is  discounted,  the  load-displacement  curves  obtained  with  larger  meshes  (h  =  13 
and  20  mm)  remain  very  consistent  with  the  X-EEM  curve  and  the  finer  mesh  results.  This  clearly 
demonstrates  the  mesh  insensitivity  of  the  new  A-EEM. 

The  mesh  insensitivity  of  the  new  A-EEM  is  further  demonstrated  by  the  fact  that  the  load- 
displacement  curves  obtained  from  unstructured  meshes  are  also  very  consistent  with  the  results 
of  the  structure  meshes  with  similar  mesh  sizes.  Two  of  the  load-displacement  curves  obtained 
from  the  unstructured  meshes,  h  =  4  and  8  mm,  as  shown  in  Eig.  11(a)  by  the  dashed  lines,  are 
basically  indistinguishable  from  the  results  of  the  respective  structured  meshes. 

Eigure  11(b)  summarizes  the  crack  trajectories  predicted  by  the  new  A-FEM  simulations 
with  the  structured  (solid  lines)  and  unstructured  (dashed  lines)  meshes.  The  crack  paths  are  all 
very  consistent  and  close  to  the  experimental  curve,  despite  they  were  obtained  with  vastly 
different  mesh  sizes  and  structures.  The  mild  difference  in  crack  paths  with  larger  meshes  is  largely 
attributed  to  the  maximum  kinking  angle  limit  of  -i-/-45°  permitted  by  our  numerical  program.  This 
is  a  standard  numerical  treatment  to  avoid  possible  crack  curve  back.  It  is  seen  that  even  with  this 
artificial  limit,  it  takes  no  more  than  2  to  3  sequential  kinks  for  the  trajectories  to  establish  a  path 
that  is  very  close  to  the  experimental  curve. 

It  should  be  mentioned  that  the  mesh  insensitivity  in  this  example  is  expected  because  all 
the  mesh  sizes  are  considerably  smaller  than  the  cohesive  zone  size  (or  fracture  process  zone  size). 
The  cohesive  zone  size,  estimated  by  =  £Tj  /  [23,  93]  is  about  =  350  mm .  This  is 

about  twice  the  beam  depth  and  more  than  15  times  the  largest  mesh  size  (20  mm).  In  this  example, 
the  mesh  size  is  limited  by  the  spacing  between  the  notch  and  the  nearest  loading  zone,  which  is 
merely  40  mm  (Eig  10a).  With  the  largest  20  mm  mesh,  there  are  only  two  columns  of  elements 
between  the  notch  and  the  nearest  loading  zone.  The  fact  that  the  new  A-EEM  simulation  can  still 
correctly  capture  the  crack  trajectory  is  rather  encouraging.  Elsewhere,  it  has  been  demonstrated 
through  other  examples  that  the  new  A-FEM  remains  fairly  accurate  with  mesh  sizes  up  to  70% 
of  the  cohesive  zone  size,  in  which  cases  the  X-EEM  exhibited  serious  (unacceptable)  mesh 
dependency  [57]. 
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Figure  11  A-FEM  simulated  (a)  load-displacement  curves  and  (b)  crack  trajectories  the  X-FEM  results  of  [66]and 
experimental  results  of  [92]  . 


Numerical  Efficiency  and  Stability  as  Compared  to  PNM-FEM 

In  this  section,  we  compare  the  numerical  accuracy,  efficiency,  and  stability  of  the  new  A- 
FEM  against  the  phantom-node  based  FEM.  The  meshes  and  boundary  conditions  are  exactly  the 
same  as  those  used  in  previous  section.  It  is  noted  that  both  the  new  A-FEM  and  the  PNM-based 
A-EEM  have  been  implemented  in  the  ABAQUS  package  (v6.11)  as  user-defined  elements  and 
they  employ  the  same  crack  tracking  algorithm  [18,21]  (the  tracking  algorithm  of  the  new  A-EEM 
was  directly  inherited  from  the  previous  PNM-based  A-EEM).  Thus  the  numerical  performance 
comparison  should  be  objective. 

For  all  simulations,  the  loading  parameters  were  set  to  be  identical  with  a  maximum 
prescribed  displacement  of  0.1  mm,  and  a  suggested  maximum  incremental  size  of  0.002  mm.  The 
minimum  increment  size  is  10'^  mm  and  largest  suggested  incremental  size  is  10'^  mm.  RIKS 
option  (arc-length  method)  was  invoked  in  attempt  to  capture  the  snap-back  behavior. 

Comparisons  of  the  predicted  load-displacement  curves  using  the  new  A-FEM  and  the 
PNM-based  A-FEM  are  given  in  Figure  12(a).  The  X-FEM  result  reported  in  [66]  is  also  included 
in  this  plot  for  comparison.  Both  the  new  A-FEM  and  the  PNM-based  A-FEM  completed  the 
simulations  without  much  numerical  difficulty.  The  new  A-FEM  results  agree  very  well  with  the 
X-FEM  result  of  Moes  and  Beiytschko  (2002).  The  PNM-based  A-FEM  also  did  a  good  job  except 
that  the  snap-backs  are  slightly  sharper  than  those  predicted  by  the  new  A-FEM  and  by  the  X- 
FEM. 

However  the  new  A-FEM  empowered  by  the  consistency-check  based  condensation 
algorithm  is  much  more  efficient  than  the  PNM-based  A-FEM.  Figure  12(b)  compares  the  CPU 
time  (in  seconds)  required  by  the  two  methods  in  completing  the  simulations.  For  an  identical 
mesh  the  CPU  time  required  by  the  new  A-FEM  is,  on  average,  merely  1/50  of  that  required  by 
the  PNM-based  A-FEM  (i.e.,  an  increase  of  numerical  efficiency  by  5,000%).  The  ~50  times 
increase  of  CPU  time  by  the  PNM-based  A-FEM  is  NOT  due  to  the  use  of  double  nodes  (recall 
that  the  total  DoFs  of  a  PNM-based  A-FEM  model  is  two  times  the  DoFs  of  the  corresponding 
model  of  the  new  A-FEM).  The  two  dashed  lines  in  Figure  12(b)  represent  the  CPU  times  required 
by  both  methods  to  complete  the  elastic  problems  with  the  same  mesh  and  the  same  load 
incremental  numbers,  but  with  crack  initiation  and  propagation  prohibited.  In  such  crack- free. 
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elastic  simulations,  the  CPU  time  required  by  the  PNM-based  A-FEM  is  about  60%  more  than  that 
required  by  the  new  A-FEM,  which  is  consistent  with  what  reported  in  [18].  Thus,  for  the  PNM- 
based  A-FEM,  the  -100  time  (10,000%)  increase  of  CPU  time  from  a  crack- free  elastic  simulation 
to  the  respective  crack  propagation  simulation  is  directly  due  to  the  additional  computational  cost 
in  treating  the  propagating  crack.  The  CPU  times  of  the  new  A-FEM  simulations  are  on  average 
about  twice  that  of  the  respective  elastic  calculations  (-100%  increase). 

A  close  examination  of  the  two  methods  shows  that  the  significant  increase  of  CPU  time 
by  the  PNM-based  A-FEM  is  associated  with:  1)  the  ever  increasing  active  DoFs  as  a  crack 
propagates,  and  2)  the  numerical  procedure  to  judicially  determine  the  crack  configuration  and  to 
properly  assign  the  ghost  and  real  nodes  to  cracked  subdomains.  These  two  elements  are  critical 
to  many  advanced  numerical  methods  that  employ  additional  external  DoFs  to  account  for 
propagating  cracks  (e.g.  X-FEM^,  PNM,  and  G-FEM).  However,  the  new  A-FEM  does  not  need 
either  of  them  (see  [21]  for  the  detailed  A-FEM  formulation).  It  is  obvious  that  the  consistency- 
check  based  algorithm  of  the  new  A-FEM  contributes  significantly  to  the  improvement  in 
numerical  efficiency  and  accuracy. 


Figure  12  (a)  Comparison  of  predicted  load-displacement  curves  using  the  PNM-based  A-FEM  and  the  new  A- 
FEM.  (b)  Comparison  of  the  CPU  time  as  a  function  of  mesh  sizes. 


6.2  Coupled  Fiber-Matrix  Interface  Debonding  and  Kinking 

We  next  demonstrate  the  capability  of  the  new  A-FEM  in  dealing  with  multiple  crack 
interaction.  We  chose  to  model  a  single  fiber/matrix  domain  as  shown  in  Figure  13.  The  entire 
domain  is  1  mm  x  1  mm  and  the  fiber  diameter  is  0.5  mm  (i.e.,  fiber  volume  fraction  19.6%).  Two 
sets  of  meshes  were  used  to  simulate  the  problem:  a)  a  fine  mesh  as  shown  in  Fig  13(a)  with  2244 
elements,  and  b)  a  coarse  mesh  as  shown  in  Fig  13(b)  with  1 10  elements.  In  the  fine  mesh  (a)  the 
elements  are  arranged  to  conform  to  the  circular  fiber/matrix  interface.  The  ring  of  elements  that 
host  the  bi-material  interface  are  initially  augmented  with  a  weak  discontinuity  and  are  allowed 


^  We  also  compared  the  numerical  efficiency  with  the  X-EEM  in  ABAQUS  (v6.1 1)  and  found  that  the  present  A- 
FEM  is  2-3  orders  of  magnitude  more  efficient  (see  94.  Mohammadizadeh,  S.,  A  novel  Augmented  Finite 

Element  Method  for  Modeling  Arbitrary  Cracking  in  2D  Solids,  in  Dept  of  Mechanical  and  Aerospace  Engineering 
University  of  Miami,  PhD  thesis,  2013.). 
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for  secondary  crack  kinking  into  the  matrix  subdomain  if  the  subdomain  stresses  meet  the 
maximum  stress  criterion.  The  coarse  mesh  (b)  is  a  structured  one  and  does  not  conform  to  the 
fiber/matrix  interface  (shown  by  the  dashed  circle  in  Fig  13(b)).  That  is,  the  fiber-matrix  interface 
is  deliberately  set  to  bisect  the  structured  elements.  The  elements  traversed  by  the  interface  are 
initially  augmented  to  account  for  intra-element  material  difference,  possible  interface  debonding, 
and  crack  kinking  into  the  matrix. 

The  modulus  and  Poisson’s  ratio  for  the  fiber  and  matrix,  respectively,  are  Ef  =  40  GPa; 
Vf  =  0.33;  and  F’m  =  4  GPa;  Vm  =  0.4.  A  cohesive  model  with  triangular  traction- separation  laws 
were  used  for  interface  debonding  and  matrix  cracking.  The  cohesive  parameters  for  matrix 
cracking  and  interface  debonding  are  set  to  be  identical:  <t  =  f  =  50  MPa; 

^nc  =  ^sc  =  1  ^  10“^  mm;  and  5^^  =  J,;  =  1  x  10  mm;  Therefore,  the  fracture  toughnesses  for  mode 
I  and  mode  II  are  Tj  =  Tjj  =  0.25  N/mm.  For  the  kinking  crack  initiation  and  propagation  in  the 
matrix,  the  maximum  principal  stress  criterion  is  assumed. 

Symmetric,  displacement-controlled  loading  is  imposed  on  both  the  left  and  right  edges, 
while  the  top  and  bottom  edges  are  stress  free.  For  both  cases,  the  prescribed  maximum  nominal 
strain  is  2%  and  the  maximum  nominal  strain  increment  is  0.05%.  Under  this  loading  scheme,  the 
CPU  time  needed  to  complete  the  simulations  is:  140  seconds  (58  increments)  for  the  fine- 
conforming  mesh,  and  20  seconds  (58  increments)  for  the  coarse-nonconforming  mesh. 

The  maximum  strain  distributions  at  four  loading  stages  for  the  two  cases  are  shown  in 
Figure  13  (a)  and  (b).  The  predicted  curves  of  nominal  stress  versus  nominal  strain  are  shown  in 
Figure  13(c).  In  both  cases,  the  nominal  stress-strain  curves  can  be  roughly  categorized  into  three 
stages.  The  first  stage  is  the  elastic  stage  (si  <  0.5%),  wherein  the  interface  remains  well  bonded 
and  the  response  is  linear,  which  can  be  seen  from  Fig  13(al)  and  (bl)^.  At  a  nominal  strain  of 
0.5%,  debonding  cracks  start  to  initiate  at  (R,  0)  and  (-R,  0),  and  they  propagates  both  upwards 
and  downwards  along  the  interface.  These  interface  cracks  propagates  in  a  symmetric  fashion 
regarding  to  both  horizontal  and  vertical  directions  (Fig  13  a2  &  b2).  As  the  debonding  length 
increases,  the  composite  stiffness  is  gradually  reduced,  but  the  stress  continues  to  rise  due  to  further 
stretching  in  matrix.  This  stage  is  eventually  followed  by  a  kinking  initiation.  For  the  fine  mesh, 
kinking  occurs  at  a  nominal  strain  of  -0.84%  and  for  the  coarse  mesh,  it  happens  at  a  strain  of 
1.2%  (Figure  13  a3  &  b3).  Shortly  after  the  kinking  initiation,  the  kinking  cracks  propagate  quickly 
towards  the  top  and  bottom  specimen  boundaries,  completely  separating  the  entire  domain  (Figure 
13  a4  &  b4).  This  process  is  characterized  by  a  sudden  drop  of  load  in  the  load-displacement 
curves. 


^  The  apparent  strange  strain  distributions  in  Fig  11(b)  are  due  to  ABAQUS  graphics  output.  Since  ABAQUS  does 
not  support  graphic  depicting  of  UELs,  we  had  to  superimpose  a  layer  of  ABAQUS  CPS4  elements  with  near  zero 
stiffness  to  display  the  strain  distribution.  In  this  way,  any  displacements  are  interpolated  continuously  within  an 
element. 
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(bl)  (b2)  (b3) 


(a)  (b) 

Figure  13  (a)  Crack  development,  including  the  initial  fiber-matrix  interface  (dashed  line),  debonding  crack,  and 
final  kinking  cracks  in  matrix,  are  indicated  in  the  top  and  bottom  contour  plots  for  the  two  meshes.  (b)The  new 
A-FEM  simulated  nominal  stress  versus  displacement  for  a  single  fiber/matrix  domain  under  uniaxial  tension. 


It  is  seen  that,  despite  the  vast  difference  in  mesh  size  and  interface  conformity,  the 
predicted  stress-strain  curves  remain  very  similar  to  each  other.  The  differences  in  predicted 
composite  elastic  modulus  and  energy  dissipation  are  well  within  5%  (The  predicted  modulus 
values,  5.2  GPa  for  mesh  (a)  or  5.37  GPa  for  mesh  (b),  are  moderately  greater  than  the  lower- 
bound  estimation  of  4.86  GPa  based  on  the  rule-of-mixture,  which  is  a  quick  proof  that  the  current 
predictions  are  reasonable).  Finally,  the  10%  and  14%  difference  in  predicted  composite  strengths 
and  ultimate  strains,  respectively,  by  the  two  meshes  are  at  least  acceptable  given  the  large 
difference  in  mesh  size  and  interface  conformity. 

It  is  emphasized  here  that  in  the  above  simulations,  the  tightly  coupled  cracking 
development  processes,  especially  the  formation  of  the  kinking  crack  from  the  active  interface 
debonding  crack  tip  zone,  is  purely  determined  by  the  program  from  the  computed  local  stresses 
in  the  matrix  subdomain  according  to  the  prescribed  maximum  principal  stress  criterion  (for 
initiation).  While  this  complex,  coupled  fracture  processes  can  be  captured  with  relative  easily  by 
the  present  A-FEM  with  subdomain  (secondary)  cracking  capability,  few  other  existing  methods 
can  truly  handle  this  with  explicit  consideration  of  the  interactive  stresses  at  the  two  active  fracture 
process  zones. 

6.3  Simulations  of  Complex  Damage  Evolution  in  a  3D  textile  composite 

In  this  example,  the  3D  A-FEM  is  used  to  simulate  the  multiple  crack  initiation  in  a  C/SiC 
angle  interlock  weave  shown  in  Eig.  14(a),  which  is  a  3D  reconstruction  from  a  high-resolution 
pCT  measurement  reported  in  [95].  A  virtual  specimen  was  generated  with  stochastic  tow 
geometry  using  Monte  Carlo  methods  in  which  the  geometrical  model  resulting  from  application 
of  the  Markov  chain  algorithm  was  converted  into  a  mesh  of  solid  computational  elements  suited 
to  the  A-FEM  formulation  as  shown  in  Eigure  14(b)  [10,  96].  The  calibrated  tow  elasticity  is 
numerically  similar  to  the  elasticity  deduced  for  a  similar  carbon/SiC  material,  reported  in  [97]. 
Nonlinear  fracture  laws  were  calibrated  from  the  same  tests  by  matching  the  predicted  matrix 
cracking  stress  with  the  experiment. 

The  matrix  containing  the  crack  is  present  in  this  image,  but  the  matrix  elsewhere  has  been 
stripped  away  in  rendering  the  image  to  reveal  the  underlying  tow  architecture.  All  such  warp  tow 
segments  are  indicated  by  the  vertical  dashed  lines  in  Eig  14(c);  they  correlate  consistently  with 
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crack  arrest.  (The  right  image  in  Fig  14(e)  actually  refers  to  a  material  in  whieh  the  weft  tow 
spacing  was  larger  than  in  the  left  and  center  images,  but  is  otherwise  of  identical  architecture;  it 
has  been  re-scaled  to  facilitate  comparisons.  Importantly,  the  relation  of  crack  position  to  weft 
tow  position  is  the  same  in  all  three  images.) 

Key  features  of  the  eraek  pattern,  ineluding  the  approximate  periodieity  of  the  mieroeraeks, 
which  form  above  near-surface  weft  tows,  and  the  alignment  of  microcracks  in  diagonal  bands 
(one  of  which  is  outlined  by  the  ellipse  in  Fig  14c),  are  well  replicated  by  the  virtual  tests.  The 
loeations  of  the  eraeks  relative  to  underlying  weft  tows  are  also  well  predieted  by  the  3D  A-FEM 
simulation.  One  partieular  eharaeteristie  of  the  mieroeraeking  demonstrates  the  potential  of  the 
A-FEM  for  virtual  tests  because  of  its  capability  in  predicting  details.  The  microcracks  seen  in 
experiments  form  above  weft  tows  and  propagate  both  into  the  composite  and  across  the  surface 
of  the  eomposite,  traeking  the  weft  tow.  However,  the  propagation  across  the  surfaee  is  limited: 
when  the  mieroeraek  approaehes  the  region  where  the  weft  tow  dips  beneath  a  warp  tow, 
propagation  is  arrested.  Thus  the  mieroeraeks  remain  of  finite  length  (Fig.  14c).  This  behavior  is 
replicated  well  by  the  3D  A-FEM  simulations.  The  images  of  Eigs.  14(c,  d)  focus  on  a  single 
mieroeraek  in  the  matrix  layer.  The  mieroeraek  ean  be  seen  to  have  followed  a  weft  tow  and  then 
arrested  as  it  approaehed  the  loeation  where  the  weft  tow  is  submerged  under  a  warp  tow. 


Figure  14.  (a)  3D  visualization  of  |iCT  data  provides  further  detail  of  the  observation  in  (a)  that  matrix  cracks 
that  have  formed  above  a  weft  tow  segment  arrest  when  they  approach  the  region  where  a  warp  tow  that  crosses 
over  the  weft  tow.  (b)  virtual  specimen  generated  from  the  pCT  data  with  top  graph  shows  the  two  architecture 
and  the  bottom  graph  shows  the  mesh  for  the  entire  model,  (cl)  textile  preform  coated  with  thin  layer  of  matrix 
(C/SiC  angle  interlock).  (c2)  observations  of  microcracks  by  digital  image  correlation  in  composite  formed  by 
full  infiltration  of  SiC  matrix  into  preform  in  (a).  (c3)  predicted  crack  pattern,  (d)  Detailed  rendering  of  a 
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simulation  shows  one  microcrack  above  a  weft  tow  correctly  predicted  to  arrest  upon  approaching  the  underlying 
warp  tow. 

The  phenomenon  of  crack  arrest  is  correlated  with  lowering  of  stress  in  the  superficial  matrix 
by  straightening  of  the  warp  tow.  (Earlier  evidence  for  tow  straightening  in  this  class  of 
composites  was  reported  in  [98].)  Both  this  stress  heterogeneity  and  the  critical  stress  for 
microcracking  depend  on  geometrical  details  such  as  the  thickness  of  the  superficial  matrix  layer 
and  the  orientations  of  tows.  They  will  therefore  be  affected  by  stochastic  variance  in  the  tow 
architecture.  By  simulating  the  evolution  of  discrete  crack  systems  with  high  fidelity,  the  virtual 
tests  can  directly  relate  variance  in  microstructure  to  variance  in  failure  due  to  microcracking  (a 
potentially  fatal  event  in  hot  structures). 

6.4  Progressive  Failure  Prediction  of  Composite  Laminates 

6.4.1  Unidirectional  Open-Hole  Tension  (OHT)  and  Compression  (OHC)  Specimens 

In  this  case  study  the  progressive  damage  evolution  in  the  unidirectional  open  hole  tension 
specimen  reported  in  Prabhakar  and  Waas  (2013)  [99]  were  simulated.  The  UD  composite  is  a 
proprietary  material  similar  to  IM7/8552.  Here  the  properties  reported  in  [99]  are  summarized  in 
Table  1,  supplemented  with  other  parameters  from  various  data  sources  from  literature.  Specimens 
with  fiber  orientation  of  90°,  45°  and  0°  were  simulated  with  both  tension  and  compression  load. 
Whenever  available,  the  A-FEM  results  were  compared  with  either  experimental  or  simulated 
results  in  literature. 


Table  1.  Composite  properties  of  the  OHT  and  OHC  specimens 


IM7/8552 

Elastic 

Constants 

^IT  ^IC 

(GPa) 

^2T  ^2C 

(GPa) 

^12 

V 

23 

(GPa) 

136 

6.7 

0.33 

0.016 

3.2 

Strength 

Xt 

X 

C 

Yt 

Y 

C 

Sl2 

(MPa) 

2960 

1680 

60 

288 

90 

Toughness 

r 

i  ft 

r 

^  FC 

r 

^  MT 

r 

^  MC 

Tn 

(N/mm) 

14.8 

25.2 

0.67 

2.47 

1.67 

90°  Tension  and  Compression 

The  90°  OH  specimen  as  shown  in  Figure  15(a)  is  50.8  mm  along  the  loading  direction  and 
38.1  mm  wide.  The  thickness  is  0.635  mm  and  the  center  hole  diameter  is  6.35  mm.  To  check  the 
mesh  sensitivity  of  the  current  method,  the  specimen  was  discretized  into  two  sets  of  unbiased, 
unstructured  meshes,  one  with  1298  elements  and  the  other  with  2093  elements. 

The  A-FEM  simulations  indicated  that  for  this  case,  crack  initiation  occurred  at  the  maximum 
stress  concentration  sites,  i.e.,  the  top-center  and  bottom-center  of  the  hole  edge.  Once  initiated. 
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both  cracks  propagated  rapidly  towards  the  speeimen  edge  leading  to  a  complete  two-part 
separation.  The  simulated  tensile  reaction  force  vs.  applied  displacement  curves  using  the  two 
meshes  are  very  consistent  and  they  are  very  close  to  the  CDFE  results  reported  in  [99].  Note 
though,  the  CDFE  results  were  obtained  with  biased  meshes,  i.e.,  assigning  much  refined  elements 
in  regions  where  the  crack  initiation  and  extension  are  anticipated.  The  predieted  peak  load 
corresponds  to  a  OHT  strength  of  25  MPa,  much  smaller  than  the  un-notched  transverse  tensile 
strength  of  Yt  =  60  MPa  (a  maximum  stress  concentration  factor  of  2.4  along  the  hole  edge). 

Under  compression  (Fig.  15b),  the  A-FEM  predicted  that  two  compressive  cracks  would  be 
initiated  at  the  same  locations  of  the  tension  cracks,  but  at  a  much  larger  (compressive)  load,  due 
to  a  mueh  larger  compressive  transverse  strength  of  Yc  =  288  MPa.  The  predicted  OHC  strength, 
deduced  from  the  compressive  load-displacement  curves  in  Figure  15(b),  is  130  MPa. 


(a) _  _ W 


(c)  (d) 

Figure  15  (a)  OHT  specimen  and  loading;  (b)  OHC  specimen  and  loading;  (c)  A-EFM  predicted  OHT-90°  load- 
displacement  curves  as  compared  to  the  CDFE  results;  (d)  A-FEM  predicted  OHC-90°  load-displacement 
curves. 


45°  Tension  and  Compression 

The  OH-45°  specimen  as  shown  in  Figure  16(a)  and  (b)  is  75.2  mm  along  the  loading 
direction  and  38.1  mm  wide.  The  thiekness  is  0.635  mm  and  the  center  hole  diameter  is  6.35  mm. 
The  specimen  was  discretized  into  two  sets  of  unbiased,  unstructured  meshes,  one  with  2441 
elements  and  the  other  with  4299  elements. 
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(a) 


—  AFE-2441  elements 
-  AFE-4299  elements 


0.5  1  1.5 

Displacement  (mm) 


(C)  (d) 

Figure  16  (a)  OHT-45°  specimen  and  loading  (2441  ele  mesh);  (b)  OHC-45'’  specimen  and  loading  (2441  ele  mesh); 

(c)  A-EFM  predicted  OHT-45°  load-displacement  curves  as  compared  to  the  CDFE  results;  (d)  A-EEM 
predicted  OFIC-45*’  load-displacement  curves. 


The  A-FEM  predicts  that  under  tension,  the  OHT-45°  was  failed  by  initiation  of  two  cracks 
near  the  top-center  and  bottom-center  of  the  hole  edge  respectively,  which  rapidly  propagated 
along  fiber  direction  (i.e.  45°)  all  through  the  specimen  width  (Fig  16a).  With  the  help  of  inertia- 
based  stabilizing  method,  this  strongly  unstable  crack  initiation  and  propagation  until  two-part 
separation  was  simulated  without  any  numerical  difficulty.  The  simulated  load-displacement 
curves  for  the  two  meshes  are  almost  identical  to  each  other  and  they  are  sufficiently  close  to  the 
experimental  and  CDFE  results  reported  in  [99]  (Fig  16c)  Furthermore,  the  predicted  crack  paths 
are  very  consistent  with  the  experimental  observation  and  the  predicted  OHT-45°  strength  of  45 
MPa  is  very  close  to  the  experimental  value  of  42  MPa  [99] . 

For  the  OHC-45°  test,  the  A-FEM  predicts  a  similar  sudden  failure  process.  However,  the 
two  cracks  are  initiated  at  different  locations,  instead  of  top-center  and  bottom-center,  the  two 
cracks  were  initiated  at  the  left-middle  and  right-middle  of  the  hole  edge  (Fig  16b)  (These  crack 
initiation  locations  were  verified  by  independent  stress  analyses  using  ABAQUS  CPS4  element 
according  to  Sun’s  criterion).  The  simulated  load-displacement  curves  for  the  two  meshes  are 
consistent  with  each  other.  The  ultimate  (compressive)  load  of  4700  N  is  much  higher  than  the 
tensile  failure  load  (1200  N)  due  to  the  larger  Yc  value.  The  predicted  OHC-45°  strength  is  196 
MPa. 


0°  Tension  and  Compression 

The  OH-0°  specimen  as  shown  in  Figure  17  is  101.6  mm  along  the  loading  direction  and  38.1 
mm  wide.  The  thickness  is  0.635  mm  and  the  center  hole  diameter  is  6.35  mm.  The  specimen  was 
discretized  into  three  sets  of  unbiased,  unstructured  meshes  with  6706  elements,  12891  elements. 
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and  26290  elements,  respeetively.  The  eharaeteristie  element  sizes  are  1.0  mm,  0.5  mm,  and  0.25 
mm,  in  the  same  order. 

In  this  partieular  ease,  since  the  experimentally  measured  load-displacement  curve  as 
reported  in  [99]  is  plagued  by  early  grip  slipping  (which  caused  premature  cracking  at  a  much 
reduced  peak  load).  We  seek  alternative  validation  from  a  cohesive  zone  model  as  shown  in  Figure 
17(b).  In  this  model,  the  two  experimentally  observed  crack  paths  were  explicitly  modelled  with 
cohesive  elements;  while  the  other  domains  are  discretized  into  standard  CPS4  elements  available 
in  ABAQUS  (V6.14). 


(b)  CPS4  +  CZM  24009  element  model 


Figure  17  OH  specimen  geometry  and  loading,  (a)  A-FEM  mesh  without  initially  designated  cracks  (12981  elements), 
and  (b)  CZM  mesh  (24009  elements)  with  two  designated  crack  paths  in  accordance  with  experimental 
observation  of  [99]. 


The  simulated  load-displacement  curves  for  the  three  different  A-FEM  meshes  and  the  CZM 
model  are  summarized  in  Figure  18(a).  In  all  cases,  the  curves  are  pretty  linear  up  to  a  load  level 
of  11000  N.  Beyond  this  point,  all  curves  exhibit  some  degree  of  deviation  from  linearity,  which 
can  be  seem  more  clearly  from  the  magnified  view  in  Figure  18(b).  Although  the  overall  load- 
displacement  responses  from  the  four  different  models  are  fairly  consistent  to  each  other,  each 
curve  exhibits  somewhat  different  local  kinks  (i.e.  sudden  load  drops)  at  different  loading  stages. 

Each  of  these  mild  load-drops  corresponds  to  the  stick-slip  behavior  of  crack  propagation, 
which  are  very  sensitive  to  mesh  size  and  local  mesh  structure  along  the  crack  path.  As  an  example, 
the  crack  initiation  and  propagation  in  the  12891 -element  A-EE  model  and  the  24009  CZM  model 
were  compared  in  Eigure  19.  Eigure  19  (al)  shows  the  contour  plot  of  8max  when  the  load -point 
displacement  is  A  =  0.355  mm,  just  before  any  significant  crack  extension  and  corresponding  to 
point  ©  in  Eigure  19(b).  The  magnitude  and  distribution  of  Smax  in  CZM  model  and  AEEM  model 
are  nearly  identical.  Note  that  at  this  point  in  the  A-EEM  model,  four  small  cohesive  cracks,  one 
for  each  quadrant  roughly  at  an  angle  of  15°  with  respect  to  the  vertical  center  line,  have  already 


40 


been  initiated  but  have  not  become  traction-free  cracks  yet.  The  predicted  initiation  locations  of 
these  cracks  are  consistent  with  experimental  observations  of  [99],  as  well  as  with  a  detailed  pCT 
characterization  of  similar  specimens  under  similar  loading  conditions  [86]. 


Figure  18  Load-displacement  relation  predicted  using  A-FEM  as  compared  with  data  from  CZM  analysis  and 
experiments. 

Shortly  after  this  point,  at  A  =  0.365  mm  (point  (D  in  Figure  19b),  one  of  the  cracks  in  the 
CZM  model  suddenly  propagated  and  arrested  at  a  distance  of  about  25  mm  to  the  left  of  the  open 
hole,  as  indicated  in  the  left  contour  plot  of  Figure  19(a2).  This  resulted  the  first  load  drop  in  the 
CZM  load-displacement  curve  in  Figure  19(b).  At  this  point  in  the  A-FEM  model  (right  plot  of 
Figure  19(a2)),  the  four  small  cracks  extended  by  a  small  amount  (~  1.0  mm),  which  caused  a 
small  reduction  of  stiffness  in  the  A-FEM  load-displacement  curve  but  no  visible  load-drop. 

Upon  further  displacement  increase  to  A  =  0.385  mm  (point  @  in  Figure  19b),  the  upper 
right  crack  in  the  CZM  model  extended  suddenly  to  about  25  mm  and  the  left  bottom  crack  also 
extended  by  a  further  amount  of  5  mm,  as  indicated  in  the  left  contour  plot  of  Figure  19(a3).  This 
resulted  a  second  load-drop  in  the  CZM  load-displacement  curve  (Fig  19b).  Between  point  (D  and 
@  in  the  A-FEM  model,  two  of  the  four  initiated  cracks,  the  upper  right  one  and  the  bottom  left 
one,  extended  almost  simultaneously  (at  A  =  0.380  mm)  by  an  amount  of  25  mm,  which  resulted 
a  significant  load-drop  in  the  A-FEM  load-displacement  curve.  The  other  two  initiated  cracks 
stopped  further  development  from  this  point  on. 

As  the  load-point  displacement  further  increased  to  A  =  0.475  mm  (point  @  in  Figure  19(b)), 
the  two  cracks  in  the  CZM  model  propagated  steadily  towards  to  the  loading  edge.  However,  in 
the  A-FEM  model,  each  of  the  two  cracks  propagated  in  a  stick=slip  fashion,  causing  several  kinks 
in  the  A-FEM  load-displacement  curve  as  seen  between  point  @  and  @  in  Figure  19(b). 
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Figure  19  (a)  comparison  of  crack  development  in  the  24009-element  CZM  model  (left)  and  in  the  12891 -element  A- 
FEM  model(right)  at  four  different  displacement  increment;  (b)  the  respective  4  instances  at  the  simulated 
load-displacement  curves. 

When  subjected  to  longitudinal  compression,  the  crack  development  is  much  different  due 
to  the  change  of  stress  state  surrounding  the  hole.  In  particular,  the  tensile  hoop  stresses  at  the  left- 
center  and  right-center  of  the  hole  edge,  coupled  with  the  low  Yt  strength  values  (60  MPa),  lead  to 
the  development  of  two  tensile  dominant  cracks  much  earlier  than  the  shear  dominant  splitting 
cracks  (Figure  20a).  The  tensile  cracks  partially  alleviate  the  stress  concentration  around  the  hole 
hence  the  development  of  the  shear  dominant  splitting  cracks  are  much  delayed. 
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(4)  A  =  0.6438  mm 


Figure  20  Strain  contour  plots  of  a  0°-OHC  lamina  under  four  load  levels  (26290  element  mesh)  illustrating  the  multiple 
crack  initiation  and  propagation  in  this  specimen,  (a)  two  tensile  crack  initiation  at  applied  displacement  of 
A  =  0.32  mm,  (b)  appreciable  propagation  of  two  additional  corner  cracks  at  A  =  0.358  mm,  (c)  much 
extended  two  corner  cracks  suppressing  the  propagation  of  the  first  two  tensile  cracks  at  A  =  0.367  mm,  and 
(d)  final  stage  of  simulation  at  A  =  0.821  mm. 
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6.4.2  Quasi-isotropic  Open-Hole  Tension  (OHT)  and  Compression  (OHC)  Specimens 

In  this  case  study,  we  apply  the  nonlinear  orthotropic  A-FEM  to  analyze  the  progressive 
damage  evolution  in  multi-directional  laminated  composites,  which  typically  exhibit  arbitrary 
multiple  cracking  events  co-evolving  from  earlier  loading  stage  up  to  final,  catastrophic  failure 
[21].  Explicitly  damage  accumulation  descriptions  of  various  damage  modes  typically  observed  in 
composites,  including  matrix  cracking,  fiber  rupture  in  tension/kinking  in  compression,  and 
delamination  are  all  included  as  described  in  section  5. 

The  quasi-isotropic  open  hole  specimen  with  a  layup  of  (45/90/-45/0)s  as  illustrated  in  Eigure 
21  is  of  dimension  150  mm  x  38  mm  x  1.0  mm  and  the  center  hole  diameter  is  6.35  mm.  Each  ply 
is  modeled  with  6707  A-EE  elements  and  all  the  interlaminar  interfaces  are  modeled  with  3D 
cohesive  zone  elements  developed  previously  in  the  Pi’s  group. 


Cohesive  interface  elements 


Figure  21  Quasi-isotropic  open  hole  specimen  with  L  =  150  mm  and  W  =  38  mm. 

The  crack  initiation  is  based  on  Kashin’s  criterion.  Upon  cohesive  crack  initiation  in  an 
element,  the  element  is  augmented  automatically  and  a  cohesive  damage  accumulation  law 
according  to  the  crack  type  is  introduced  to  account  for  the  post-crack  fatigue  life.  Eor  the  material 
studied  in  this  paper  (IM7/8552),  the  cohesive  damage  accumulation  laws  for  all  the  crack  types 
have  been  previously  calibrated  and  reported  by  [100,  101]. 

Eigure  22  gives  a  validated  example  of  the  A-EEM  prediction  to  an  open-hole  tension-tension 
specimen  reported  in  [102].  The  specimen  is  exposed  to  a  cyclic  loaded  with  a  maximum  load 
level  of  50%  of  the  failure  load  and  a  load  ratio  0.1  for  30,000  cycles,  followed  by  a  monotonic 
load  until  complete  failure.  Eigure  22(a)  shows  the  strain  contours  in  the  plies  as  well  as  the  fatigue 
cracks  generated  automatically  by  the  A-EEM.  The  A-EEM  predicted  cracks  and  crack  types  are 
consistent  with  experimental  observation.  The  predicted  residual  strength  of  5 15  MPa  is  very  close 
to  the  experimentally  reported  mean  value  of  544  MPa.  The  deamination  at  the  interfaces  as  shown 
in  these  figures  is  also  consistent  with  experimental  observations  [102].  This  and  other  validated 
examples  indicate  that  the  A-EEM  based  fatigue  analysis  tool  can  meet  the  challenge  of  predicting 
fatigue  life  of  composites  with  explicit  consideration  of  multiple  discrete  damage  events. 
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Figure  22  (a)  A-FEM  predicted  individual  cracking  in  each  ply  (white  bands)  and  delamination  at  the  three  interfaces 
(white/red  areas)  at  the  peak  load;  (b)  much  extended  ply  cracking  and  delamination  immediately  after  peak 
load 

6.5  Coupled  BM-AFEM  for  Large  Textile  Composites 

In  this  section,  we  demonstrate  the  useful  application  of  the  developed  breakable  ID  tow 
elements  (section  5.4)  in  modeling  the  progressive  damage  in  large  textile  composites.  Shown  in 
Figure  23(a)  is  an  integrally  woven  CMC  panel  of  dimension  of  5 1  mm  x  5 1  mm  x  2.4  mm.  The 
panel  consists  of  103  weft  tows  and  99  warp  tows,  part  of  which  can  be  seen  from  the  magnified 
view.  The  material  properties  of  this  panel  has  been  calibrated  from  a  leveraging  AFRL  program 
and  will  not  be  repeated  here. 

The  fiber  tow  grid  and  elements  are  generated  using  the  geometric  model  generators  as 
described  in  section  5.4.  The  panel  is  discretized  into  8  layers  of  3D  solid  elements  (8  x  10,000  = 
80,000)  which  are  connected  by  7  layers  of  3D  CZM  interface  elements  (7  x  10,000  =  70,000).  In 
addition,  50,000  tow  elements  are  embedded  in  the  3D  model  with  some  tow  elements  passing 
through  the  various  interfaces. 
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Figure  23  (a)  BM-AFEM  model  for  the  integrally  woven  CMC  panel  with  embedded  tow  elements  showing  by  the 
magnified  view;  (b)  BM-AFEM  simulated  surface  stress  distribution  very  consistent  with  experimental  measurement. 
The  magnified  view  of  the  upper-right  corner  show  the  cracked  interface  bridged  by  tows,  (c)  BM-AEEM  simulated 
stress  (normalized)  vs.  strain  curves  with  initial  flaws  of  various  sizes. 

Initial  flaws  of  size  1%,  9%,  16%,  and  25%  with  respect  to  the  entire  area  (5 1  mm  x  51  mm) 
is  introduced  into  the  center  interface  to  investigate  the  effects  of  such  flaws  on  the  structural 
integrity,  especially  the  through-the-thickness  strength.  To  mimic  the  test,  the  bottom  surface  is 
fixed  and  the  entire  top  surface  is  under  a  displacement  controlled  loading. 
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The  BM-AFEM  simulated  deformation  and  surface  stress  distributions  are  shown  in  Figure 
23(b).  The  pseudo  periodic  stress  distribution  with  oscillating  tensile  and  compressive  stress, 
despite  the  overall  tensile  loading,  is  very  typical  of  textile  composites  [97,  103].  It  is  mainly 
caused  by  the  reinforcing  tow  architecture.  The  magnified  view  clearly  shows  the  tow  bridging 
across  the  debonded  interface,  which  improves  the  damage  tolerance  of  the  relatively  brittle 
CMCs. 

The  predicted  stress-strain  curves  (normalized)  are  summarized  in  Figure  23(c).  All  stress- 
strain  curves  exhibit  an  initial  spike  followed  by  a  prolonged  plateau  with  about  10%  of  the  peak 
strength.  The  initial  spike  with  strain  less  than  0.01%  is  due  to  the  very  brittle  nature  of  the  matrix 
(with  critical  energy  release  rate  on  the  order  of  1-10  J/m^):  the  interface  debonding  occurs  at  very 
small  strain  levels  (<  0.01%).  However,  the  fiber  tows  are  able  to  bridge  the  debonded  interface 
and  provide  further  damage  tolerance.  It  also  shows  that,  for  this  particular  loading  condition, 
effective  bridging  is  not  fully  realized  until  the  nominal  strain  reaches  about  3-5%.  Other 
simulations  show  that  fiber  bridging  is  much  more  effective  when  the  fracture  process  is  shear 
dominant. 

It  is  rather  encouraging  that,  although  the  fiber  tows  are  modelled  as  the  simplest  possible 
ID  elements  (which  are  embedded  in  the  3D  domain  with  multiple-point-constraint),  the 
simulation  results  are  nevertheless  very  consistent  with  physical  reality. 


7.  CONCLUDING  REMARKS 

In  this  study,  a  new  augmented  finite  element  method  (A-FEM)  has  been  formulated  based 
on  the  principle  of  virtual  work.  One  of  the  major  advantages  of  this  method  is  that  it  can  accurately 
account  for  arbitrary  intra-element  cracks  and  their  interaction  without  the  need  for  additional 
external  DoFs  as  in  X-FEM  or  extra  nodes  as  in  PNM.  The  new  A-FEM  formulation  does  not  need 
to  assume  deformation  modes  for  elemental  displacement  enrichment  as  in  the  embedded 
discontinuity  method  either.  Instead,  it  introduces  internal  node-pairs  with  normal  displacements 
as  internal  nodal  DoFs,  which  are  eventually  condensed  at  elemental  level.  Thus  the  crack 
displacements  are  natural  outcomes  of  the  elemental  equilibrium  consideration.  The  advantage  is 
multi-fold.  First,  it  enables  a  unified  treatment  of  both  weak  and  strong  discontinuities.  Second, 
for  strong  discontinuities  with  piece- wise  linear  constitutive  relations,  the  new  A-FEM  formulation 
ensures  an  exact  solution  (in  the  piece-wise  linear  limit)  to  the  condensed  elemental  equilibrium, 
which  greatly  improves  the  numerical  accuracy,  stability,  and  efficiency  by  2-3  orders  of 
magnitude.  Third,  the  unique  A-FEM  formulation  makes  it  straightforward  to  repeat  the 
augmentation  procedure  within  an  element  to  include  multiple  intra-element  cracks,  which  is  very 
powerful  in  dealing  with  multiple  crack  interaction  problems  because  such  augmented  elements 
do  not  need  extra  external  DoFs  or  nodes. 

The  elemental  condensation  was  done  through  a  novel  consistency-check  procedure.  It  has 
been  shown  that  it  provides  analytical  solution  to  local  equilibrium  for  piece-wise  linearized 
cohesive  laws.  Rather  than  assuming  trial  crack  displacements  and  iterating  for  elemental 
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equilibrium  such  as  in  the  Newton-Raphson  method,  it  starts  with  trial  cohesive  segments  in  the 
cohesive  laws  and  finds  the  analytical  solution  (in  piece-wise  linear  sense)  to  the  elemental 
equilibrium  through  a  simple  consistency  check.  For  piece-wise  linear  cohesive  laws  with  only  a 
small  number  of  possible  stiffness  segments,  this  algorithm  is  very  efficient. 

Through  the  numerical  examples,  it  has  been  demonstrated  that  the  new  A-FEM  method, 
empowered  by  the  new  elemental  condensation  algorithm,  achieved  very  2-3  orders  of  magnitude 
improvement  in  numerical  accuracy,  efficiency,  and  stability.  The  A-FEM’s  capabilities  in  high- 
fidelity  simulations  of  interactive  cohesive  cracks  in  heterogeneous  solids  have  also  demonstrated 
through  the  benchmarking  examples  in  this  report  (more  substantial  results  are  published  in  journal 
articles  acknowledged  to  this  project  as  listed  in  appendix  A).  Therefore,  it  is  concluded  that  all 
the  proposed  research  goals  as  listed  in  the  beginning  of  this  report  have  been  achieved. 

Eurthermore,  we  have  gone  beyond  the  proposed  goals  of  this  study  to  expand  the 
capability  and  to  apply  the  A-EEM  to  benefit  several  other  projects,  through  a  wide  collaboration 
with  national  and  international  university,  US  government  laboratories,  and  US  industrial 
companies  (See  details  in  Appendix  C).  These  ongoing  collaborations  have  generated  sufficient 
interest  elsewhere  to  be  continuing  beyond  the  life  of  the  present  program. 

On  the  educational  front,  the  project  has  benefitted  three  (4)  PhD  students,  two  (2)  MS 
students,  and  two  (2)  post-doctorate  researchers.  Most  of  them  have  graduated  and  joined  the  US 
workforce  as  researchers  in  US  research  institutions  and  industrial  companies. 
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Model  (BM),  two  strategies  for  textile  CMCs  have  been  tested:  1)  coupled  BM-AFEM  that 
can  explicitly  account  for  the  matrix  cracking,  and  the  temperature-  and  environment- 
dependent  tow  failure;  2)  coupled  BM-CDM  that  uses  calibrated  continuum  damage 
mechanics  model  (CDM)  to  account  for  the  material  degradation  due  to  numerous  initial 
flaws  and  matrix  microcracking.  The  first  strategy  is  far  more  computationally  costly  but  can 
provide  high-fidelity  material  failure  information  (suitable  for  microscopic  or  unit  cell  level 
analysis  with  adequately  characterized  initial  material  configuration).  The  second  strategy  is 
less  accurate  and  requires  experimental  calibration  from  physical  testing  of  real  a  material, 
but  it  does  provide  a  quick  and  effective  means  for  rapid  analyses  of  large  CMC  structures 
(see  section  6.5  for  example).  This  is  a  key  contribution  to  the  research  community  in  ultra- 
high  temperature  (UHT)  materials  and  structures. 

2.  DARPA  Materials  Development  Platform  (MDP),  PI:  Olivere  Sudre,  Teledyne  Scientific 
Company. 

The  thermal-mechanical  A-EEM  (with  steady-state  and  transient  temperature  DoEs)  have 
enabled  us  to  team  up  with  Teledyne  (led  by  Dr.  Olivier  Sudre)  to  test  its  capability  in  virtual 
testing  of  new  materials  of  interest  to  DARPA.  The  coupled  BM-CDM  strategy  serves  as  an 
essential  platform  in  virtual  development  and  assessment  of  high-temperature  composite 
structures  such  as  airfoils.  The  program  is  expected  to  transfer  to  Teledyne  and  DARPA  by 
December  2017. 

3.  Micromechanical  Failure  Analysis  of  Aligned  Low-Aspect  Short  Fiber  Composite 
Materials  for  TFF,  PI:  Ryan  Karkkainen,  University  of  Miami.  DARPA  Award  No. 
HROOI 17258 15. 
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In  this  project,  the  3D  A-FEM  is  being  used  to  assess  the  mechanical  properties  of  short- 
fiber  reinforced  polymers.  The  fact  that  the  A-FEM  allows  for  explicit  representation  of  all 
failure  modes  at  microscopic  level  (fiber  breakage,  matrix  cracking,  and  fiber/matrix 
interface  debonding)  in  a  material  domain  of  significance  volume  may  offer  a  new  way  of 
predicting  such  composite  properties  without  homogenization,  which  typically  requires 
calibration  parameters  from  real  materials. 

4.  Composite  Progressive  Damage  Analyses.  The  Boeing  Company.  PI:  Qingda  Yang. 

We  have  been  supported  by  Boeing  Structural  Damage  Technology  Group  since  2014  to 
further  develop  the  A-FEM  for  progressive  fatigue  damage  analysis  with  explicit 
consideration  of  multiple  discrete  crack  formation  and  propagation.  Under  this  support,  the 
2D  A-FEM  has  been  extended  greatly  to  include  many  material  characteristics  crucial  to 
composite  failure,  including  material  orthotropy,  shear  nonlinearity,  temperature  dependent 
material  and  fracture  properties.  More  critically,  the  2D  A-FEM  now  has  the  capability  of 
initiating  multiple  types  of  cracks  observed  in  different  plies  in  a  composite  laminate  (matrix 
cracking  in  tension,  shear,  and  compression,  fiber  rupture  in  tension  and  kinking  in 
compression,  all  coupled  with  possible  delamination  as  shown  in  section  6.4.2),  under  either 
static  or  cyclic  loadings.  The  2D  fatigue  A-FEM  has  been  successfully  transfer  to  Boeing 
and  is  currently  being  subjected  to  their  internal  evaluation  (with  their  proprietary  material 
properties).  Boeing  is  currently  supporting  us  to  develop  a  3D  brick  element  which  can  be 
more  convenient  for  laminated  composites. 

5.  Continued  Collaboration  with  Army  Research  Laboratories  (ARL) 

Since  the  beginning  of  this  project,  we  have  been  in  collaboration  with  researchers  in  ARL 
(Dr.  C.F.  Yen)  to  work  on  problems  of  mutual  interests.  During  the  project  period  two  of  our 
students,  Derek  Schesser  (then  PhD  candidate,  2014)  and  Scott  Floerke  (then  MS  student, 
2017)  worked  with  Dr.  Yen  as  summer  interns.  Derek  helped  to  transfer  the  earlier  version 
of  the  2D  A-FEM  into  the  ARL  computing  platform  and  worked  under  Dr.  Yen’s  guide  on 
the  development  and  instrumentation  of  a  novel  high-speed  torsional  test  facility  at  ARL. 
Scott  Eloerke  worked  on  the  dynamic  fracture  of  aluminum  panel,  using  both  our  dynamic 
A-EEM  and  the  commercial  software  LS-Dyna  and  ABAQUS.  This  collaboration  resulted  a 
MS  thesis  and  a  journal  paper  is  in  preparation.  We  are  ready  to  transfer  the  much  more 
matured  version  of  the  A-EEM  upon  request. 

6.  Collaboration  with  National/International  Universities 

We  have  been  in  close  collaboration  with  Dr.  Spearing’s  group  in  Southampton  University 
(UK,  Dr.  Mark  Spearing)  on  application,  validation,  and  verification  of  the  A-EEM 
predictions  with  in-situ  micro-CT  measurements  of  composites.  Dr.  Spearing’s  group  has 
access  to  the  state-of-art  micro-CT  measurement  facilities  and  has  a  wealth  of  in-depth 
material  toughening  data  and  their  correlation  with  the  in-situ  observation  of  the  microscopic 
damage  development.  This  collaboration  has  so  far  resulted  in  two  high-profile  journal 
publications  and  this  collaboration  is  continuing. 

We  are  also  in  collaborating  with  Professors  in  Demark  Technical  University  (DTU,  Profs. 
Bent  Sorensen,  Brian  Legarth)  on  coupling  A-EEM  simulation  with  their  state-of-art 
capability  of  in-situ  SEM  for  microscopic  measurement  of  fiber/matrix  bonding  properties. 
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The  major  research  findings  have  been  published  in  the  ASME  Journal  of  Applied 
Mechanics. 

In  collaboration  with  the  Department  of  Engineering  Mechanics  at  Zhejiang  University, 
we  have  studied  the  failure  mechanisms  in  syntactic  foams  and  the  cyclic  thermal-mechanical 
response  of  fiber/matrix  interfaces,  which  resulted  two  joint  publications. 

We  are  also  in  collaboration  with  a  visiting  student  from  the  Northwestern  Polytechnical 
University  of  China  (Mr.  Junsan  Hu)  to  investigate  the  joining  technologies  for  composite 
materials.  In  particular,  we  attempt  to  gain  understandings  through  experimental  testing  on 
the  local  damages  surrounding  those  pin-holes  in  typical  bearing  joints  due  to  mechanical 
insertion  and  in-service  attrition.  We  have  co-authored  two  journal  publications. 

Dr.  Brian  Cox,  our  subcontractor  at  Arachne  Consulting,  has  been  collaborating 
extensively  with  Dr.  Yang  and  with  third  parties  in  three  projects  that  were  undertaken  to 
broaden  the  scope  of  applications  of  AEEM  in  problems  of  interest  to  the  Army  and  the 
materials  development  community.  We  have  been  in  collaboration  Teledyne  to  Generating 
code  for  textile  reinforced  components;  with  Dr.  Jerry  Quek  of  the  Institute  for  High 
Performance  Computing  (IHPC),  Singapore,  to  explore  methods  of  generating  defect 
populations  at  the  fiber  scale  in  a  nominally  unidirectional  ply  by  extremely  fast 
computational  methods  based  on  graph  theory;  with  Dr.  X.  Boyang  in  the  Delft  University 
(Netherland)  to  study  the  fundamental  question  of  how  bifurcation  of  a  crack,  or  the 
coalescence  of  two  cracks,  occurs  when  the  mother  crack  and  the  mother/daughter  pair  are 
governed  by  a  nonlinear  spatially-distributed  fracture  zones. 

All  of  these  collaboration  projects  have  generated  sufficient  interest  elsewhere  to  be  continuing 
beyond  the  life  of  the  present  program. 
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